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I.  INTRODUCTION  AND  MOTIVATION 

A  well-documented  module  ETBFCT  for  solving  generalized  continuity 
equations  was  presented  in  NRL  Memorandum  Report  3237,  dated  March  1976. 

The  module  PRBFCT  was  also  included  to  treat  the  case  of  periodic  boundary 
conditions.  In  these  modules,  the  cell  centers  are  specified  while  the 
cell  boundaries  are  located  midway  between  the  centers.  In  August  1980, 
JPBFCT,  in  which  the  cell  boundaries  are  tracked,  was  documented. 

The  above  modules  are  based  on  the  Flux-Corrected  Transport  (FCT) 
technique  introduced  first  by  Boris  and  Book.1'  FCT,  instead  of  adhering 
to  an  asymptotic  ordering,  requires  positivity,  a  physical  and  mathematical 
property  of  continuity  equations.  To  assure  positivity,  the  convective 
stage  includes  or  is  supplemented  by  a  large  diffusive  flux  of  zeroth 
order  (in  s  =  •  Consequently,  an  antidif f usive  or  corrective  step  has 

to  follow.  The  two  stages  together  are  able  to  tre-  .  steep  gradients 
without  generating  dispersive  ripples.  Antidiffusion  being  a  physically 
(and  numerically)  unstable  process,  the  corrective  flux  is  limited  according 
to  a  criterion  which  may  be  stated,  "The  antidif fusion  stage  should  generate 
no  new  maxima  or  minima  in  the  solution,  nor  should  it  accentuate  already- 
existing  extrema. " 

FCT  was  shown  to  be  applicable  to  any  finite  difference  transport 

2 

scheme  and  able  to  improve  it.  Phoenical  FCT,  a  refinement  which  minimizes 
residual  diffusive  errors,  was  introduced.  Clipping  and  terracing,  two 
nonlinear  processes  resulting  from  the  flux  limiter  were  discussed. 

Finally,  splitting  techniques  were  recommended  to  extend  FCT  to  multi¬ 


dimensions. 


2 

The  most  detailed  error  analysis  of  FCT  algorithms  was  performed  in 
Ref.  3.  Low-residual-diffusion  and  lew-phase-error  algorithms  were  derived. 

An  optimal  algorithm,  Fourier  FCT,  was  introduced. 

The  requirements  for  positivity  of  a  general  three-point  scheme  and 
the  antidiffusion  flux  for  a  minimum  residual  diffusion  were  derived  in  Ref.  4. 

Zalesak^  provided  a  general  mathematical  interpretation  of  the  anti¬ 
diffusion  flux  as  the  difference  between  a  high-order  transport  scheme  and 
a  low-order  one.  He  also  described  a  generalized  fully  multidimensional 
flux  limiter  guaranteeing  that  the  antidiffusion  fluxes  on  all  sides  of  the 
control  volume,  acting  in  concert,  do  not  create  any  ripples.  It  was  shown  that 
by  proper  selection  of  the  flux  limiter  parameters  the  clipping  and 
terracing  phenomena  can  be  reduced. 

The  goal  of  the  present  work  is  to  extend  JPBFCT  to  a  fully  two- 
dimensional  algorithm,  without  time  splitting,  and  incorporate  the  Zalesak 
flux  limiter  while  still  keeping  the  implementation  of  the  convective, 


diffusion  and  antidiffusion  processes  as  physical  fluxes. 


II.  FOURIER  ANALYSIS;  DEFINITIONS 


3p 

at 


A  generalized  conservation  equation  can  be  written  in  the  form: 
+  U.Vp  =  pV.U  +  s (x, t, 0/ • • • ) / 


(1) 


where  u  is  the  velocity  vector,  P  is  the  generalized  density  or  the  trans¬ 
ported  quantity  whose  positivity  is  to  be  conserved,  and  s  is  a  source  term 
including  all  the  remaining  terms,  i.e. ,  gradients,  divergences,  body  forces, 
etc. 


In  the  analysis  we  assume  s  =  0  and  u  =  constant.  We  shall  i  *t  with 
the  one-dimensional  case.  Eq.  (1)  reduces  to 

It  +Uo  It  =  °'  (2) 

whose  analytic  solution  is 


p (x,t)  =  p  (x  -  u  t,0) , 
o 


(3) 


a  rightward-propagating  wave  with  velocity  u  -  Let  us  Fourier  analyze 


j  (x,t)  in  space,  assuming  periodic  boundary  conditions.  Assuming  an  initial 
distribution  of  density  p(x,0)  =  F(x), 


F  ( 


r — * 

x)  =  —  +  )  (a^  cos  kx  +  sin  kx) , 


(4a) 


k=l 


2:7  k 


where  k  is  assumed  to  be  normalized,  i.e.,  k  replaces  — ; —  .  In  complex 
form 

oo 

_ ,  .  V""k  *  ikx 

F0°  "L  Vs 

k=~® 

where  i  = 

quantity  pkis  related  to  these  by 

3 


(4b) 


yj~- 1.  From  the  reality  of  p(x,0),  the  and  b^  are  real.  The 


4 


f3k  "  lbk 


for  k  ■>  0; 


a,  +  10. 

,  X  K  _ 

=  < - r -  tor  .<  <  0; 


(5) 


0  for  k  =  0. 


Notice  that  we  could  have  started  the  summation  in  Eq.  (4a)  from  k  =  0  since 

sin  0=0  and  cos  0=1.  The  zeroth  order  term  would  then  be  a  .  The  form 

o 

(4a)  is  preferred,  however,  since  it  is  compatible  with  the  symmetric  formul 

tion  of  Eqs.  (4b)  and  (5).  Then  ;  is  given  by 

K 


x ) e  ax . 


(6) 


•ron  Eq.  (3) ,  the  density  prefile  at  time  t  is  given  by 


r— 'la:  .  ik(x  -  u  t)  ikx 

:u't)=L  :ke  °  =L  vt)e  - 


(7) 


vnere 


Vfc)  =  V 


-iku  t 


(8 ) 


shewing  that  each  harmonic  independently  advances  uniformly  in  phase  without 

changing  its  magnitude  (see  Fig.  1). 

Suppose  p  is  known  at  all  times  only  as  a  set  of  N  +  1  quantities 

:  ;  on  discrete  grid  points  with  separation  ox  "  jr;  x  =  j  ox  (j  =  0,1,..., 

J  N  3 

M 

N-l),  since  .  =  c  .  We  can  have  only  —  +  1  different  narmonics.  Namelv, 

O  N  2 

wave  numbers  (0 , 1 ,  .  .  .  ,N/2)  and  wavelengths  (°°  =  L/0,  y,  of,  .  .  .  , 
respectively,  where. we  note  that  the  shortest  wave  length  :  s  26 x.  Let 

Ao  r-' N/2 

f  vx)  =  —  +  J  (A^  cos  kx  +  sin  kx)  ,  (9a) 

k=l 


*  - 


or 


5 


f  (x) 


-E 


N/2 


.  ikx 


(9b) 


k=-N/2 


2ti(N/2) 

(see  Fig.  2).  We  notice  in  Eq.  (9a)  that  at  k  =  N/2,  sin  kx .  =  - - -  ;Px 

3  b 


sin  tt j  =  0;  hence  only  cos  kx  is  needed  at  k  =  N/2.  There  are  then 


N  coefficients,  A^  (k  =  0,...,N/2)  and  B^  (k  =  1,2,..., N/2-  1),  which  can 


be  determined  using  f(x^)  =  p  (j  =  0,...,N-1),  where  superscript  0  denotes 


2  it  -N  -  2+  n 

time  t  =  0.  Similarly,  since  for  all  j,  exp  [i  — -  (~)jox]  =  exp  (i—  (— )j:x 

Ld  2.  Li  L 


Eq.  (9b)  is  rewritten  as 


f(x)  =^2 

k- -N/2+1 


N/2  1  ikx 


pke 


(9c) 


Again,  we  get  N  coefficients  5^(k  =  -N/2+1, ... ,0, ... ,N/2) .  The  relation 


between  the  p^  and  A^,  are  given  by  equations  similar  to  Eq.  (5). 
Formally , 


=  Ar 

N  L-J 


N-l 


3=o 


o  ikx  . 
P.e  3 


(10) 


Eq.  (3)  predicts  the  density  at  time  t  as 


?  (x,t)  = 

k=-N/2+l 


N/2  3  eik(x  -  UQt)  ^  N/2  2  ikx 

k  =  )  k 


-E 

k=-N/2+l 


where  p,  (t)  =  P.e  Since  we  are  only  concerned  with  P(x.,t), 

K  K  J 


substituting  x  =  x^  =  iox,  we  get 


5:v°  -E 

k=-N/2+ 1 


N/2  -  ikj5x 
Pk(t)e 


(ID 


r 


T 


If  the  time  is  also  discretized,  let  tn  =  nft,  o* 
then 


(x  ,tn) ,3?  =  i,  (tn)  , 
:  k  x 


‘J-E 


N/2  ~n  ikjfx 
“k6 


(12) 


k=-N/2+l 


where 


=  * 


P,  e 


-iku  nit 


(13) 


If  we  space-discretize  only,  after  we  Fourier  analyze  the  initial  density 

profile  p°,  i.e.,  after  gettina  the  ,5,  in  Eq.  (9c),  the  problem  is  reduced 
1  K 

iicx 

to  that  of  propagation  of  the  complex  harmonics  e  (k  =  0,...,N/2).  In 
a  nonlinear  problem,  each  harmonic  can  couple  into  components  of  the  other 
harmonics.  In  the  linear  problem  of  Eq.  (2) ,  however,  each  harmonic  propagates 
independently  (this  is  also  true  if  u  =  u(t)).  Since  the  number  of  spatial 
points  does  not  change,  we  can  always  express  the  density  at  any  time  as 
a  Fourier  expansion  of  the  form  Eq.  (9c).  In  a  nonlinear  problem  p^tt) 


is  a  function  of  (5 


-N/2+1 


,o  , .  . .  ,3  .  ,J  at  time  t  =  0.  But  in  the  linear 

o  N/^. 


problem  o^tt)  is  only  a  function  of  as  is  obvious  from  Eq.  (13) . 
the  time  is  also  discretized,  we  can  then  define  a  transfer  function 


A  (k)  = 


„n+l 

fk _ 

„n 

o. 

K 


(14) 


which  is  independent  of  n  if  u  =  uq  as  is  obvious  from  Eq.  (13)  (analytic 
solution) ,  yielding 


-iku  ft 
A  (k)  =  e  o 


Eq.  (12)  may  be  rewritten  then  as 

n  ikjfx 


1-E 


N/2  3k[A(k)l“  e' 


(15) 


(16) 


k=-N/2+l 


*-r  ryr.i  ’’B* |  . 


i 

j 


7 


u  6t 

Denoting  the  constant  ^  by  £  and  the  dimensionless  wave  number  k-'x  by  £, 
“  i  B  ^  1 

A ( S)  =  e  .  The  amplification  is  | A  ( 3 ) |  =  1  and  the  phase  shift  is  -Be. 

2r  N 

Notice  that  the  smallest  £  =  0  and  largest  8  =  —  —  f x  =  t.  For  a  finite- 

1j  2. 

difference  scheme  applied  to  the  linear  problem,  each  harmonic  propagates 
independently.  Consequently,  a  method  equivalent  to  Fourier-analyzing 


n+l 


and  pn(j  =  0,...,N-1)  and  evaluating  A.  from  Eq.  (14)  is  to  study  the 
13  K 

_  ,  ,  .no  ikj£x  ,  o  . 

propagation  of  only  one  harmonic  by  assuming  =  o  e  ,  where  p  is 

constant.  Then 


A  (k)  = 


n+l 

J _ 

n 

P1 


(17) 


3y  writing  A(k)  as 


A  =  |A|ei’3, 


(18) 


we  define  the  amplitude  (or  diffusion)  error  and  the  relative  phase  error 


as 


a  =  At  -  1 


(19a) 


and 


D  _  (-9)  ~  Sc  _  ~(9/8)  . 

Sc  c 


(19b) 


respectively.  We  define  a  scheme  as  stable  if  A  <  1  (see  Fig.  3). 


Example : 


Assuming  u  =  const  =  u  ,  the  original  explicit  SHASTA  Algorithm  can  be 


written  as 


TD  n  e ,  n  n  .  ^  ,1  ,  e  .  ,  n  _  n  n  . 

j  j  2  j+1  i-1  v8  2  vkj+l  pj-l' 


n+l  TD  1 ,  TD  _  TD  TD  , 
°i  ■  °i  -  -  20  j  +  °j-i> 


(20) 


8 


in  which  we  identify  -  j(p”+1  -  as  the  net  transportive  flux,  denoted 

by  a  superscript  T,  and  +  cr  )  (o'1,,  -  2pn  +  o11  )  as  the  net  diffusive 

*  8  Z  J  +  1  J  J  —  1 

TD 

flux,  denoted  by  a  superscript  D,  from  which  we  have  the  notation  p 
TD 

Expressing  p  as  a  three-point  formula  we  can  write 


TD  ,1  £  e,  n  -,i  ,  e  n  r  e  e  u 

°j  =  [8  +  I  "  2Ipj+l  +  [1  -  2{8  +  2  )lpj  +  [8  2  +  I^j-l 


e,  n 


Each  of  the  quantities  in  square  brackets  is  >_  0  for  j  e  j  <_ -r ,  assuring  the 
TD  n 

positivity  of  o ,  if  a .  >_  0.  The  positivity  requirement  will  be  discussed 

,  .  ,  . ,  .no  ikjcx 

later  in  detail.  Assuming  p .  =  p  e  , 


TD  o  ikj<5x  e.  o  ik(j+l)6x  o  ik(j-l)<5x 
=  o  e  -  -(p  e  -  o  e  ) 


+  (8  +  T  >  (p 


o  ik(j+l)5x  o  ikj5x  o  ik(j-l)5x, 
e  -  2p  e  +  p  e  )  , 


giving 


TD  .  n  .  e  ,  iS  -iS,  .  ,1  ^  e*\  .  i£  .  -i3 

p .  / p .  =  1  -  —  (e  -e  )+(-+— )(e  -2+e  ). 

3  3  -  a  z 

i3  -i3  . ,  .  - 

e  -  e  ib  -ip 

Denoting  the  operator - — - -  =  i  sin  g  by  t  and  e  -2  +  e 

TD  n  1  c 4 

=  2  (cos  3-1)  by  d,  we  have  p .  =  (1  -  et  +  vd)  c .  where  v  =  —  +  — 

3  384 

_  i  4,  pn+1  =  (1  -  Et  +  vd)  on  -  ud(l  -  et  +  vd)pn,  whence 
83  3  3 


Then  i 


A  = 


n+1 

P  . 

"1 

- d. - 

n 

D  . 

J 


=  (1  -  et  +  vd) (1  -  ud) 


(21) 


We  notice  that  at  e  =  0,  A  /  1.  In  fact,  A  =  (1  +  -rd)  (1  -  — d)  ,  a 

o  o 

deficiency  that  led  to  the  introduction  of  a  phoenical  algorithm  in  Ref.  2, 
in  which  the  antidiffusion  operates  on  a  transported  density  which  is  free 
from  any  zeroth-order  diffusion.  Phoenical  SHASTA  is  written  as 


£»*>+  t. 


9 


T  n  e ,  n  n  e  n 

pj  =  pj  -  2"  P j+1  '  pj-l}  +  2  (pj  +  i 


-  2pn  +  ,n  _> 

3  3-1 


TD  T  1  ,  n  _  n  n  , 

°J  '  °j  *  3toj+l  '  2CJ  *  “j-l” 


n+1  TD  1 ,  T  ,  T  T  . 

>5  *  Cj  -  atpj+i  -  p  p3.i» 


(22) 


thus  yielding 


A  =  (1  -  et  +  Xe“d) (1  -  ud)  +  vd. 


(23) 


where  \  =  v  =  y  =  satisfying  A  =  1  at  e  =  0. 

2  0 


The  importance  of  phoenicity  lies  in  the  fact  that  the  total  diffu¬ 
sion  through  a  surface  is  proportional  to  the  time  of  diffusion  and  there¬ 
fore  should  vanish  as  6t  -+  0,  i.e.,  c  0 . 

Later,  in  Ref.  6,  ETBFCT  and  JPBFCT,  based  on  the  scheme 


T  n  e ,  n  n 

Pj  -  Pj  -  2(pj+1  -  Pj_1) 


TD  T  .  n  _  n  n  , 

Pj  =  Pj  +  V(pj+1  "  2pj  +  °j-l} 


n+1  TD  ,  T  _  T  T  , 

Pj  =  Pj  -  y(Pj+1  -  2P j  +  Pj.,), 


(24) 


were  introduced,  yielding 

A  =  (1  -  £t) (1  -  yd)  +  vd,  (25) 

1  2  ,  £2 

where  v  =  —  +  —  and  y  =  —  -  —  .  Notice  that  the  zeroth  order  term  is  the 
6  3  6  6 

2 

same  in  both  v  and  v,  thus  yielding  a  residual  diffusion  0(e  ),  which 
vanishes  as  6 1  -*•  0. 


*"*  ■  w- 


•Vf 


III.  AMPLITUDE  AND  PHASE  ANALYSIS 


If  in  Eq.  (18)  A  is  expressed  as  A  =  A  +  iA  ,  where  R  stands  for 

R  X 

real  and  I  for  imaginary,  then 

i  1 2  2  2 

|A|  =  Ar  +  Ax,  (31a) 

0  =  tan  1 (A  /A  )  .  (31b) 

X  R 

Equations  (31)  yield  numerical  values  of  [a(  and  0  for  a  given  8-  These 
should  be  expanded,  however,  in  a  power  series  in  6  and  plugged  into  Eas. 

(19)  to  get  an  estimate  of  the  order  of  a  given  scheme.  Expanding  Eas.  (31) 
in  power  series  is  a  huge  task.  Instead,  we  use  a  scheme  based  on  successive 
differentiation,  as  follows: 


PHASE  ERRORS 


As  seen  from  Eqs.  (21) ,  (23)  ,  and  (25)  ,  three-point  schemes  can  be 
expressed  in  terms  of  a  transport  operator  t  =  i  sin  2  and  a  diffusion 
operator  d  =  2  (cos  8-  1).  In  other  words,  A  =  A(t,d)  where  t  =  t(8)  and 
d  =  d(S).  Taking  the  logarithm  of  Eq.  (18),  we  obtain  log  A  =  log  | A j  +  i0, 
yielding 

9  *  Imflog  A] .  (33) 

Expanding  9  in  a  power  series  of  B,  near  S  =  0,  we  have 


9  =  9  + 

o 

where  (  ) '  i 
from  Eq.  (33) 


p*  I 
o  1 ! 


d  (  ) 
dS 


and  the  subscript  0  denotes  the  value  at  8 


0.  Since, 


10 


11 


1109  a,1s-°h 

all  we  need  are  the  derivatives  of  (log  A)  with  respect  to  6,  at  3  =  0. 

First,  by  direct  differentiation  we  get  (log  A)'  =  A'/A, 

2 

(log  A)  "  =  A”/A  -  (A'/A)  and  so  on.  Noticing  that  the  "consistency"  of 
any  scheme  requires  A(B  *  0)  =  1,  we  can  write 

(log  A)q  =  Aq;  (35a) 

(log  A)o  =  Aq  -  Aq2;  (35b) 

(log  a/"-  <"  -  3aV  ♦  2A^3,  ,3501 

»V  t  v  ff«»  rr  2  *2M  '4 

(log  A)  =  A  -  4A  A  -  3A  +  12A  A  -  6A  ;  (35d) 

O  O  OO  O  OO  o  '  / 

V  v  I  I  V  Iim  1 2  "•  noi 

(log  A)  =  A  -  5A  A  -  10A  A  +  20A  A  +  30A  A 

OO  OO  OO  oo  oo 

’3  "  ’5 

-  60Aq  Aq  +  24Aq  .  (35e) 

Next  denoting  — by  (  )  and  ■—  (  )  by  (  )^(  we  get  by  direct  differentia- 
tion  A  =  t  A  +  d  A  ,  A  =  t  A  +  d  A  +  t  A  +2td  A  +  d  2A,  and 

so  on.  Confining  our  scope  to  schemes  of  first  degree  in  t  (composite 
transport  excluded)  and  of  second  degree  in  d,  we  have 
. tt  .  tdd  ddd 

A  =  0,  A  =  constant,  and  A  =0.  (36) 

We  obtain  then 


Ao  =  1? 


'  '  t  '  d 

A  =  t  A  +  d  A  ; 
o  o  o 


A,,  =  t"At  +  d"Ad  +  2t,d,Atd+  d’2Add; 
o  oo  oo  ooo  oo 


1  "1 
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•  ••  "■  t  d  "  '  *  "  td  '  ”  dd  1  '2  tdd 

A  =  t  A  +  d  A  +  3  (t  d  +  t  a  )A  +3ddA  +3td  A  ; 

o  oo  oo  oo  ooo  ooo  ooo 

i  V  i  V  f  i  V  n  till  Hi*  i  in  j-j 

A  =t  A^+d  A  +  (4t  d  +  6t  d  +  4t  d  )A 

o  oo  o  oo  oo  ooo 

i  in  n  t  i  i  ii  ii  i2  tdd 

+  (4d  d  +  3d  ) A  +  (12t  d  d  +  6t  d  )A 

OO  O  Q  OOO  oo  o 


V  V  t  v  d  'V  I  Ml  I.  It  in 

A  =  t  A  +  d  A  r  (5t  d  +  lQt  d  +  lot  d 

OO  OO  OO  OO  OO 

I  I  V  td  '  '  V  "  dd  II  1  II 

+  St  d  )A +  (5d  d  +  lOfl  d  )A  +  (30t  d  d 

OO  O  OO  OOO  OOO 


i  "2  ii  mi  mi  i2  tdd 

+  5t  d  +  20t  d  d  +  lOt  d  )A  . 
oo  ooa  ooo 


Going  back  to  the  definition  of  t  and  d 


»  l«  Ml  IV  V 

t  =0,t  =i,t  =0,t  =  -i ,  t  =0,  and  t 

o  o  o  o  o 


»  «i  in  iV  v 

d  =  0,d  =  0,d  =-2,d  =  0,d  =  2,  and  d 

o  o  o  o  o  o 


Substituting  in  Eqs .  (37),  wa  get 


A  =  1  ; 
o 


1  t 

A  =  lA  ; 
o  o 


A  =  -2A  ; 
o  o 


A  =  -i(Av  +  6A  ) ; 
o  o  o 


'V  d  dd 

A  =  2  (A  +  6A  )  ; 
o  o  c 


AV  =  i{At  +  30Atd  +  60Atda) . 
ooo  o 


Finally,  with  Eqs.  (30),  Eqs.  (35)  yield 


»  ■ -f  *»  - 
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(log 

A)o  " 

0; 

(log 

A)  = 

iAt; 

'40a) 

o 

o 

(log 

1 1 

A)  * 

-2Ad  +  (aV; 

(40b) 

o 

o  o 

in 

t  d  td 

t  3, 

(log 

A)o  = 

-iA  (1  -  6A  )  -  i  6A 
o  o  o 

+  2(A)  ); 
o 

(40c) 

•V 

dd  d  d 

t  2 ,  t ,  td 

(log 

A) 

■  12A  +  2(1  -  6A  )  [A 

-  2(AU)  ]  -  6A  [4AtQ 

o 

o  o  o 

o  o  o 

+  (Afc)3]; 
o 

(40d) 

v 

tdd  t  r ,  d 

,  dd  d  2 

(log 

A)  = 
o 

60iA  +  iA  [1-30A  - 

o  o  o 

60A  +  120  (A  ) 

o  o 

+  20i(At)3[l  -  6Adj  +  30iAtd[l  -  4Ad  +  4(Afc)2] 
o  o  o  o  o 

+  24i(Afc)5.  (40e) 

o 


Eqs.  (40)  invoke  the  fact  that  only  the  odd  derivatives  of  log  A 
are  imaginary.  Therefore,  with  the  use  of  Eqs.  (33)  and  (34) ,  we  get 

■  ii< 


9 


(log  A) 
_ o 

i 


(1°g  A)0  ^ 

i  3 ! 


(41) 


Example: 

Let  us  phase-analyze  the  scheme  described  by  Eqs.  (22),  i.e.,  the 
transfer  function  of  Eq.  (23) 


U  =  ~  and  X  = 

8  2 


A  =  (1  -  et  +  Xe  d) (1  -  ud)  +  vd 
where  v 

First  we  notice  that  it  is  phoenical:  A  =  1  at  e  =  0.  By  direct 

t  d  2  2 

differentiation,  A  =  -  e(l  -  ud)  and  A  =  Xe  (1  -  ud)  -  u(l  -  et  +  Xe  d)  +  v 


Atd  =  CU,  Add 


nx  2  a  7Ytdd  n 

2\z  Ur  ana  A  =  0. 


/  - 
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,  2 
AC 


~> 


t  d  2  c~ 

At  6  =  0,  t  =  0  and  d  =  0,  yielding  A  =  -  z,  and  A  =  Xc  +  (v  -  u)  =  —  ; 

_  o  o  2 


.  td 


A  *  cu  =  §■  /  Ad<^=  -  2Xe“  u  =  -•§•,  and  At"'d  =  0.  Substituting  in  Eas.  (40)  , 

O  o  O  bO 


(41) ,  we  get 


9  =  -  e3  +  7(7  e  -  e3)63  + 

6  4 


Using  Eq.  (19b) ,  the  relative  phase  error  is  found  to  be 


R  =  |(J  -  £2)32  +  0( B4) , 


showing  that  the  scheme  is  second  order  in  phase. 

Alternatively,  let  us  derive  an  expression  for  v  which  renders  the 
scheme  of  Ecs.  (24)  fourth-order  in  phase.  Upon  differentiating  the  trans¬ 
fer  function  A  *  (1  -  ct) (1  -  yd)  +  vd  we  get,  when  we  substitute  3=0, 


m  t  ,  d 

A  =  -  e;  and  A  =  v  -  p ; 
o  o 


td  dd  .  ,tad  . 

A  -  cu,  A  =0,  and  A  =0, 
o  o  o 


(42) 


which  with  Eqs.  (39a)  through  (39c)  gives 


A  =  -  1£ 

o 


A  =  -  2(V  -  y) 
o 


A  =  -  i  (-e  +  6su  )  . 
o 


(43) 


Substituting  in  Eq.  (35c) ,  we  obtain 


3  2 

Im  (log  A)  =  e  (1  -  6y)  -  6e(v  -  y)  +  2e  =  e(l  -  6v  +  2 e  )  . 


,2  . 


To  reduce  the  coefficient  of  6  in  the  R  expansion  to  zero  we  require 


(log  A)  *  0,  yielding 


J  =  6  +  3 


* 


(44) 


AMPLITUDE  ANALYSIS 

Denote  the  complex  conjugate  by  a  bar  on  top: 

| A j “  =  A  A. 
dn  _  dn 

Since  -  (  )  =  — - -  ,  we  get  by  successive  differentiation  of  Eq.  (45) 

dgn  d6n 

( [  A[  2)  =  1; 

1  1  o 

,  .  2  1  ~ 

(A  )  =  A  A  +  A  A  ;  ! 

o  o  o  o  o 

y  m  *«  »  »  rr 

(|A  )  =  A  A  +  2A  A  +  A  A  ;  ( 

o  o  o  o  o  o  o 

Iff  Ml  _  II  T  i  li  ill 

(|a|‘)  =A  A+3AA+3AA  +AA  ;  ( 

1  1  o  oo  OO  OO  CO 


9  i  V  r  v iii  »  1 1  ii  tin  iV 

( I A  “ )  =  A  A  +4A  A  +6AA  +4AA  +  A  A  . 

o  oo  oo  oo  oo  oo 


Noticing  from  Eqs.  (39)  that  the  odd  derivatives  of  A  are  pure 
imaginary  while  the  even  ones  are  pure  real. 


A  =  +  A  =  1; 
o  o 


A  =  -  A  ; 
o  o 


A  =  +  A  ; 
o  o 


III  III 

A  =  -  A  ; 
o  o 


»  V  f  V 

A  =  +  A 
o  o 


a-.  . 
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Substituting  in  Eqs.  (46) ,  we  get 


C  S  a| 

2)o  =  1; 

(48a) 

( |  a) 

2  1 

)  =  0; 

(48b) 

o 

1 

A 

i  > 

•y  1 1  II 

o  2 

(|a[ 

“)  =  2  [A  + 

(-r)  ]; 

(48c) 

o  o 

1 

2 

( |  A ! 

>o  =0; 

f  It 

(48d) 

o  i  v  i  v 

A  A 

o  o  ' 1  2 

(|a| 

“)„  =  2 [A  + 

o  o 

4  (— r-  )  (— — )  +  3  (A  )  ]  , 

11  o 

(49e) 

where  we  notice  that  the  odd  derivatives  vanish.  Accordingly,  |a|“  can  be 
expanded  as 


|A|2  =  1  +  (jAl2) 


i  i  i2  ' v  S ' 

t  +  <  W  )c  |t  - 


(49) 


Example: 

Let  us  derive  an  expression  for  u  to  render  the  diffusion  error  of 
ETBFCT  fourth  order.  Substituting  Eqs.  (42)  into  Eq.  (39d) ,  we  find 


a'V  =  2  (v  -  u)  ,  (50) 

o 

,  ,  2  "  ? 

Using  Eqs.  (43)  with  (48),  we  obtain  ( ( A |  )  =  2[-2(v  -  u)  +  c  ],  which  has 

to  vanish  for  a  fourth-order  diffusion,  yielding 


v  -  u 


(51) 


Solving  Eos.  (44)  and  (51),  we  have  u  =  —  ,  whence  —  = 

6  l 


-  e,  A_  =  -  s 


o  3  ,v  2 

— : —  *>  e  (1  -  6u)  =  £  ,  and  A  =  c  .  We  can  then  write 
l  o 


i  i  2  '  V  2  3  4  2  2 

(  A|  )  =  2[e  +  4  (-£)  (e  )  +  3e  ]  =  2e  (1  -  z  )  , 

1  1  o 
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which  when  substituted  into  Eq.  (49)  gives 
2 

|  A  |  2  =  1  yq-  (1  -  e2)«4  +  0<|3b)  ,  (52) 

4 

showing  a  slight  instability  near  g  =  0  (the  coefficient  of  g  is  positive) . 

A  warning  is  in  order  at  this  point.  Although  a  positive  coefficient 
of  the  leading  term  in  the  expansion  imp.'ies  unstable  behavior,  a  negative 
one  does  not  guarantee  a  stable  scheme,  since  the  expansion  is  valid  only 
near  3=0. 


Figure  4  shows  the  amplification  |a|  versus  g.  We  notice  a  maximum 
value  of  j  A |  =  1.0018  at  3  =  53.668°  ±  0.001  for  e  =  -g.  We  can  get  rid  of 
the  potential  instability  by  using  a  slightly  different  expression  for  u, 


By  trial  and  error,  a  was  found  to  be  1.056.  The  dashed  line  in  Fig.  4 

shows  the  resulting  amplification  for  a  =  1.056.  The  maximum  value  of  | A  j 

becomes  0.999998  at  5  =  45.775°  +  0.001.  Since  the  phase  error  depends  on 

v  only,  the  resulting  scheme  is  still  fourth-order  in  phase  error.  The 

zeroth-order  antidiffusion  being  kept  at  — ,  phoenicity  is  preserved,  i.e., 

6 

2 

the  residual  diffusion  is  0 (c  ).  Later,  a  modified  algorithm  which  is 
stable  and  has  sixth-order  diffusion  and  fourth-order  phase  error  is 
described. 


“  >~r 

-■«-*  "1 _ 
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IV.  POSITIVITY  AND  ANTIDIFFUSION 

The  concept  underlying  FCT  is  "positivity."  This  means  that  the 
sign  of  the  dependent  variable  must  be  preserved  under  the  influence  of 
convection  alone.  Source  terms  can  alter  the  sign.  Positivity  is  particu¬ 
larly  important  near  steep  gradients  where  the  convective  fluxes  tend 
to  make  the  transported  quantity  undershoot  or  overshoot.  Positivity  is 
ensured  by  supplementing  the  convective  step  with  a  large  diffusive  flux 
of  zeroth  order  in  it.  For  example,  in  the  scheme  of  Eq.  (24),  consider 
the  transport  step  alone, 

T  n  s ,  n  n  , 

h  ‘  h  *  5«Vi  ‘  °j-i) ' 


applied  tc  the  discontinuities  of  Fig.  5(a)  and  (b) ,  where  £  =  +  1/2.  The 
negative  density  in  Fig.  5(a)  and  overshoot  in  Fig.  5(b)  are  obviously  major 

errors.  P.y  supplying  enough  diffusion, 

o 

TD  T  ,1  s“  .  n  _  n  n 

h  ■  (5  *  3  -  2ci  *  h-i’' 

we  see  the  negative  density  in  Fig.  5(a)  disappear,  as  does  the  overshoot 
in  Fig.  5(b).  Formally,  in  the  expression 

o™  =  tl  -  2(~  +  f  )  ]  o'?  +  i  (m  +  ~  )  -  #]  Pn  +  [  (f  +  f  )  +  |]  Jn  , 

1  6  3  ]  6  3  2  3+1  6  3  2  j-l 

the  quantities  in  square  brackets  are  all  >  0  for  je|  <_  1/2,  therefore 

TD  n 

ensuring  positivity  of  o  .  as  long  as  o  .  ^  0. 

3  3  — 

A  side  benefit  of  the  zeroth-order  term  is  more  accurate  propagation 

i.e. ,  high-order  phase  preservation.  As  seen  from  Eq.  (44),  selecting 
^  ^.2 

v  =  —  +  —  assures  a  fourth-order  phase  error. 

A  byproduct  of  this  large  added  diffusion  is  antidiffusion,  which 

is  needed  to  extract  at  least  the  zeroth  order  part.  This  leaves  a  residual 
2 

diffusion  0(e“)  near  almost  uniform  distributions.  Near  steep  gradients, 
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antidiffusion  fluxes  have  to  be  reduced  enough  to  maintain  the  positivity 
TD 

of  a  .  This  process  is  called  correction  of  fluxes,  and  gives  rise  to 
the  name  "f lux-corrected  transport."  In  the  case  of  a  discontinuity ,  the 
local  antidiffusion  flux  is  cancelled  completely.  This  trimming  means 
that  the  amplitude  no  longer  has  the  order  of  accuracy  derived  above. 

But  near  steep  gradients  the  concept  of  order  is  meaningless  anyway.  On 
the  other  hand,  Eq.  (44)  is  independent  of  u.  The  fourth  order  phase  error 
is  therefore  assured  regardless  of  the  antidiffusion  fluxes.  Specifically, 
"the  antidiffusion  stage  should  generate  no  new  maxima  or  minima  in  the 
solution,  nor  should  it  accentuate  already  existing  extrema"  (Ref.  1). 

The  first  mathematical  formulation  of  the  above  statement  was  given 
in  connection  with  explicit  SHASTA,1 


TD 

o  . 

1 


f  }  {p 


n 

j+i 


(61a) 


n+1 

0  • 


~TD 

3 


(61b) 


The  corrected  antidiffusion  flux, 


f.  ,  =  sign  A.  ,  •  max  {0,  min  [A.  ,  •  sign  A.  ,, 
D+i  *  3+i  D-i 


■i", 

8  '“j 


.  ,  I  ,  A .  „  •  sign  A  .  . ] } 

3+i  3+3/2  *  3+t 


(62) 


is  the  corrected  form  of  the  raw  flux 


f  .  =  |-  A.  . 

3  +  i  8  3+j 


1 ,  TD  TD. 
8(Pj+!  ~  Pj  5' 


(63) 


which  in  this  scheme  is  always  in  the  same  direction  as  the  gradient  in 
TD 

o  .  There  are  eight  different  possible  cases,  shown  schematically  in 
Fig.  6.  Cases  5-8  are  mirror  images  of  1-4,  respectively. 


vf 
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Equation  (62)  will  cancel  an  antidiffusion  flux  whenever  it  would 
lead  to  accentuate  a  maximum  or  a  minimum,  as  illustrated  in  Fig.  6,  and 
will  trim  it  enough  not  to  generate  a  new  maximum  or  minimum  whenever 
it  is  not  cancelled. 

Later  in  Ref.  2  the  raw  antidiffusion  fluxes  were  evaluated  using 

T  _  1  T  T 

.  in  the  raw  flux  f.  ,  =  —  (p  .  ,  -  a  . )  ,  where 
i  T+i  8  i+l  i 


s  .  n 
2 (3j+l 


n 


)  . 


The  corrected  flux  is  expressed  as 


f.  ,  =  sign  f.  •  max  {0,  min  (1,  .  •  sign  f. 

3  +  i  3  +  i  J~t  J+i 


ifj*r-Vv2  '  si,n  V*11' 


(64) 


where  we  get  sixteen  possible  cases  (twice  as  many  as  before,  depending 


whethei 


:  .  ,  is  parallel  to  2  ,  or  coposite  to  it).  In  Fig.  7  we  consider 

j+i  3+i 


only  those  cases  when  f^  +  ^  is  positive,  since  the  other  cases  are  their 
mirror  images. 

Again,  the  flux  is  cancelled  whenever  it  would  accentuate  a  maximum 

or  minimum.  But  it  is  also  cancelled  in  cases  6-9  where  it  would  net  in 

general  cause  any  problems,  an  unnecessary/  action.  This  is  due  to  the 

fact  that  f .  .  is  corrected  independently  of  f .  ,  and  f.  „ 

j-i  3+3/2 

Zalesak^  reexpressed  the  role  of  the  flux  limiter  as  "guaranteeing 
that  the  two  antidiffusion  fluxes  associated  with  each  cell,  acting  in 
concert,  should  not  create  any  ripples."  The  mathematical  formula 
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implementing  the  above  statement  is  described  for  1-D  schemes  by  the  follow¬ 
ing  steps : 

p*  5  sum  of  antidif fusive  fluxes  "into"  grid  point  j 


=  max  (0,f.  .)  -  min  (0,f.  ,) 
3-i  3  +  4 


+  _  ,  max  TD, 

Q .  =  (C  .  -  0  .  ) 

3  3  3 


min  (1,Q+/P+)  if  P+  >0 
3  3  3 


if  P .  =0. 
3 


Similarly , 


=  sum  of  antidif fusive  fluxes  "out  cf"  grid  point  j 


max  (0 , f  ^ )  -  min  (0,f_. 


(cTD  -  omin) 
j  3 


min  (1,Q./P.)  if  P.  >  0 
3  3  3 


if  P .  =  0, 
3 


,  max  ,  min  .  n+1  ^ .  ,  ... 

where  o  .  and  c  .  are  the  upper  and  lower  oounds  on  p  .  ,  respectively,  wnicn 

.j  3  3 

ensure  that  no  ripples  form  at  grid  point  j.  Defining  the  correction  ratio 


in  (R .  ,R . )  if  f  >0 
3+1  3  3  +  4 


in  (R+,R  .  )  if  f  .<  0. 

3  3+1  3+4 


we  set 


f  =  C.  ,f.  .. 
3+4  3+4  3+4 


-9 


served  after  the  antidiffusion  step  is  performed. 
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i.  a 

A  =  -  z,  A  =  (v  +  Ae  )  -  y  / 
o  o 


td  dd  2  _tdd  n 

A  =  eu,  A  =  -2Ae  y,  and  A  =  0 
o  o  o 


Then 


A  =  -  i£; 
o 


A  =  -  2  [  (v  +  Ae  )  -  u] ; 
o 


A  =  -  i  (-  e  +  6  eu]  ; 
o 


*  v  -7  2 

A  =  2[v  +  Ae“  -  u  -  12.X  e  vi  1 ; 


For  a  fourth-order  diffusion  error. 


(|a|2)o  «  2  [-2  (v  +  Ae  -  y)  +  £  ]  =  0, 


yielding 


,  2  e 

v  +  Ae  -  y  =  — 


ti  »v 

Going  back  to  the  A  ,  A  expressions,  we  cm  rewrite  them  as 


o  o 


Aq  =  -  £  , 


•v  2  2 

A  *  e  -  24  Ae  u. 
o 


For  a  fourth-order  phase  error 


♦  ll  p 

(log  A)  =  -i[-c  +  6ey  ]  -  3(-ie)(-2)(v  +  Ae 
o 


y)  +  2(-ie) 


=  ie [1  -  6 (v  +  Ae2)  +  2e2]  =  0, 


yielding 


v  +  Ae2  =  1/6  +  £2/3, 


«  -  -  *•- 
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which  gives 


u  =  1/6  -  i"/6 

Mf  »v 

We  can  then  rewrite  A  and 

o  ° 


A 

o 


.  3 

=  it  / 


Checking, 

('a|2)'V=  2te2  -  4\e2(l  -  e2)  -  4e4  +  3e4] 

1  1  o 

=  2  [1  -  4 X )  £ 2  ( 1  -  £2)  ]  . 

showing  that  we  can  make  the  scheme  sixth-order  in  diffusion  by  selecting 
X  =  1/4. 


In  summary. 


v  =  i/6  +  £2/12 ,  u  =  1/6  -  e2/6,  X  =  1/4- 


(84) 


Again,  we  have  to  check 

2 


TD 


TD  r,  _.l  e“,  n  .1  £  £,  n  rl  £  e.n 

=  [1  -  2 (-  +  -)} a  +  -  +  -  -  -] 0 . ,  +  tr  +  T  +  ~] ° ■  .• 

j  6  2  3  6  3  2  j  +  j.  6  3  2  3-1 


Each  quantity  in  square  brackets  is  0  if  |ej  <_  1/2,  yielding  >  0  if 
0^  >  0,  thus  ensuring  positivity.  Figure  8  shows  j  A  j  and  R  versus  3. 
Finally,  we  note  that  this  is  still  a  5-point  scheme. 


V-  * 
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VI.  EXTENSION  TO  HIGHER  ORDERS  IN  DIFFUSION 
AND  PHASE  ERRORS 

We  seek  a  combination  of  transport  operator  t(t  5  i  sin  8)  and 

diffusion  operator  d(d  =  2  (cos  8  -  1))  which  approaches  the  analytic  solution 

up  to  a  prescribed  order  of  8.  Since  the  transfer  function  of  the  analytic 

-iSe 


solution  is  expressed  as  A  =  e 
A  =  cos  Se  -  i  sin  8e 


(85) 


or 


A^  =  -  sin  Be 


A  =  cos  Be* 
f\ 


Now,  we  write  sin  Be  as 


sin  Be  =  sin  8  [A  +  A. (1  -  cos  8)  +  A_ (1  -  cos  3)  +...], 

o  1  2 


(86) 


(87) 


where  Aq,  A^ ,  A^  ,  * . .  are  determined  such  as  to  make  the  series  expansion 
of  both  sides  of  Eq.  (87)  agree  up  to  a  prescribed  order  of  £.  In  other 
words,  the  derivatives  of  both  sides  with  respect  to  3  at  8  =  0  have  to 
be  equal.  We  get  the  following  system  of  algebraic  equations: 


-e 


1 

1-1 


0 

3 

-15 


0 

0 

30 


0 

0 


(88) 


solved  by  "forward  substitution"  since  the  matrix  of  coefficients  is  already 
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" lef t-triangular. "  We  solve  for  Aq  first,  then  for  A,,  etc.  We  get 


A  ,  ,  A  =  iiklij  A  ,  £(4-t2)  (1-e2) 
o  A1  3  '  2  30 


(89) 


As  for  the  construction  of  the  matrix,  the  first  column  is  the 

odd  derivatives  of  sin  8,  the  second,  those  of  sin  8(1  -  cos  8),  the  third, 

2 

those  of  sin  8(1  -  cos  8)  ,...  and  so  on,  all  at  8=0.  We  notice  that 

the  even  derivatives  are  all  zero.  To  get  these,  let  <j>  =  1  -  cos  g,  and 

define  K  recursively  K.  ,  =  K.?  where  K  =  sin  8 •  If  we  have  the  deriva- 

1  +  1  x  o 

tives  of  K. ,  those  of  K  ,  will  be 
l  l+l 


K.  ,  =  K.i; 
l+l  1 


Ki+1  -  V  +  V  ; 


I I  II 


K  .  ,  =  K  .  <f>  +  2K.  $  +  K.  <p  , 
l+l  l  1  1 


(90) 


and  so  on.  Generally  if  (  )  ^  =  - — ^-j  and  (  )  ^  =  (  )  ,  we  get 


d8  8=o 


K(n)  =y-n  n  o 
l+l  L _ ’  m  i 


(m)  (n-m) 


(91) 


m=0 


i-  ,n,  n! 

where  (  ) 


,  ,  . ,  ,■  .  All  we  need  then  is  the  derivatives  of  K  I  sin  3 

m  (n-m) !m!  o 

and  Oil-  cos  3  at  8=0;  namely. 


K  =  0,  K  =  1,  K  -  0 ,  K  =  -1,...,  and 
o  o  O  o 


<j)  =  0 ,  ^  =  0 ,  $  =1,  ■)>  =0,  V  =  - 1 ,  •  •  • 


(92) 


Now,  we  write  cos  8e  as 


cos  3e 


B  +  B  (1  -  cos  8)  +  B„,(l  -  cos  8)2  +  ... 


(93) 
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where  B  ,  B, ,  B„  , . . . ,  are  determined  such  as  to  make  the  series  expansion 
oX2 

of  both  sides  of  Eq.  (93)  agree  up  to  a  prescribed  order  of  6.  In  this 
case 


K 

o 


1,  k' 

o 


I  f 


K 

o 


=  0. 


Using  Eq.  (91) ,  we  get 


(94) 


(95) 


where  we  notice  again  a  "left-triangular"  coefficients  matrix.  By  "forward 
substitution"  we  obtain 

2  -.-2  2 

B  =  1,  B  =  -e  ,  B  =  —  (1  -  e  )  ...  (96) 

o  1  2  6 

Obviously,  we  can  get  Eq.  (93)  by  differentiating  Eq.  (87)  and 
vice  versa ,  but  we  need  then  to  continue  the  expansion  one  more  term  and 
use  trigonometric  identities.  The  direct  approach  followed  is,  however, 
preferred,  since  it  enforces  3  given  form  on  the  expansion  which  is  in 
no  way  unique,  as  explained  below. 

Noticing  that  sin  6  =  t/i  and  1  -  cos  B  =  -d/2,  we  can  write  a 
sixth-order  diffusion  error,  sixth-order  phase  error  scheme,  for  example, 


1 


as 


»R  ■  1  *  T  *  -  2411  -  e‘,a' 

iAj  -  -  et(l  -  d  ♦  i-^(l  -  ^)d2]. 


If  we  step  at  A^ ,  B^,  we  get  ETBFCT,  which  has  fourth-order  diffusion 

•) 

and  phase  error.  Although  t  and  d  are  both  three-point  operators,  td~ 
is  a  seven-point  formula.  An  important  conclusion  follows:  We  need  three 
points  for  a  second-order  diffusion  and  phase  error,  five  points  for  a 
fourth-order  error,  and  so  on,  adding  two  points  at  a  time.  We  can,  however, 
get  sixth-order  diffusion  and  fourth-order  phase  accuracy  with  only  five 

,  ,  •>  7  2 

points  since  we  have  to  match  the  sum  ;A|  =  A"  +  A  up  to  a  prescribed 


order  of  3  and  not  A  and  A  separately.  Scheme  (82)  is  an  example.  Alter- 

I  R 


natively  one  can  construct  a  scheme  with  fourth-order  diffusion  and  sixth- 


order  phase  accuracy  using  only  five  points  since  we  have  to  expand  tan 


-1  AI 


not  A^.  and  AR  separately. 

Before  implementing  Eq.  (97) ,  it  is  important  to  emphasize  that 
the  expansion  is  not  unique.  For  example,  we  can  use  the  expansions 


sin  Sc  =  sin  S  [Aq  +  A^l  -  cos  8)  +  A 0  (1  -  cos  2 3)  + 


cos  Sc  =  B  +  3. (1  -  cos  0)  +  B„(l  -  cos  23)  + 
o  1  2 


We  get  then 


10  0 


-1  3  12 


1  -15  -120 


31 


By  solving  the  two  systems  (99) ,  (100) ,  we  obtain 

,  £  (9  -  e2)  (1  -  e2)  .  -e(4  -  c2)  (1  -  -;2) 

Ao  =  £'  Al  =  15  '  A2  =  60  (100) 

and 


B  =  1 ,  B .  =  ^  (4  -  £  2)  ,  B  =  \  (1  -  f-  )  .  (101) 

o  13  2  2  6 


We  notice  that  the  matrices  are  full  and  the  coefficients  (A^ , 

A  9  » . . . ) ;  (B^,  , . . . )  are  more  complex  in  form  than  the  corresponding 

coefficients  Eq.  (89),  (96).  Moreover,  they  change  if  the  expansion  is 

extended  to  higher  order.  The  operator  (1  -  cos  23)  results  from  a  five- 

point  formula;  namely,  o.  .  -  2p .  +  o.  _.  It  is  abandoned  therefore  in 

3+2  3  1-2 

favor  of  the  t.  ae-point  operator  formula  of  Eq.  (97)  since  the  latter 
requires  knowledge  of  only  one  point  outside  the  boundary. 

We  rewrite  Eq.  (97)  as 


6 


Noting  that  pn  2  =  Apn,  let  us  collect  the  terms  in  such  a  way  as  to  ensure 

n  T  n 

positivity  at  every  step.  First,  (-st)c  =  a  -  p  where 


T  n  t  ,  n  n 

h  '  cd  •  2l5j+r  h-i- 


2  2  2  2 

n+1  n  ,e  n  ,,1  -t  ..  ,e  ,  n  ,1  -s 

:  =  P  +  (-  d)  a  -  [  ( — - —  a)  (—  d)  ]p  +  ;1-  ( — - —  d) 

2  b  4  b 

.  2  2 

.1  -  £  ,,  rlM  T  n, 

+  ( — 7 -  d)  [-(1  -  --  ) d]  }  (p  -  p  )  , 

o  5  4 


whence 


i  ->  2  ~> 

T  .  3  -  e“  ,  n  ,1  -  e  ,,  ,  T  .  ,  E  n 


5  o  +  [(-r  +  — - )dio‘  -(— 

2  o  6 


d)  (p  +  d)  p  ] 


2  2 

+  {  l1  d)  tpd  -  f  )d] }  (pT  -  tn) 

•3  5  4 

From  our  earlier  experience,  the  combination 

T  ,£2  1  -  t2,  ,  n 

0  +  (-  +  — g - )  ao 

is  known  to  be  positive  for  |e|  <_l/2.  The  remaining  terms  are 
regarded  as  antidiffusion.  The  following  scheme  is 
recommended: 


T  n  £ ,  n  n  . 

Pj  =  °j  '  I(pj  +  1  -  °j-l); 


TA  T  1  . .  £  T  T  T 

P  .  -  p  .  -  ~(1  -  7  )  P  .  ^  -  2p  .  +  0  .  .  ; 

1  J  $  4  1+1  1  1-1 


(104b) 


TAD  TA  1,.  2.  ,  n  n  n  , 

).  =  P .  +  -(1  +  £  ) (p  -  2d .  +  D .  ; 

3  3  5  1+1  J  1-1 


(104c) 


TD  T  1  £(,n  _  n  n  , 

o,  =  o  .  +(t  +  7  )  (P._, .  -  2p  .  +  p.  )  ; 

1  1  o  3  j+1  3  2-1 


(104d) 


n+1  TD  ,1  -  £  ,  TAD  „  TAD  TAD  * 

°j  =  °j  "  )(pj+i  -  2oj  +  cj-i}  ' 


(104e) 


^  -t 


where  the  asterisk  of  Eq.  (104e)  means  that  the  antidif fusive  fluxes  in 


TAD 

this  step  are  to  be  corrected.  It  is  worth  noticing  that  if  o  is  taken 
T 

as  p  ,  we  obtain  a  fourth-order  diffusion,  fourth-order  phase  algorithm. 

If 

2 

TAD  _  T  E  n 
p  =  p  +  —  do  , 


we  get  a  sixth-order  diffusion. 


TAD  _  T  1  , ,  2,  ,  n 

5  =  p  +  j(l  +  s  )do 


fourth-order  phase,  and  finally, 


yields  a  fourth-order  diffusion,  sixth-order  phase  error  scheme.  The  ampli¬ 
tude  and  phase  error  versus  3  are  shown  in  Fig.  9. 
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VII.  PHYSICAL  ASPECTS 


The  conservation  of  mass,  momentum,  and  energy  applied  to  a  system 
are  expressed  as 


0 


V£  (t) 


_d_ 

dt 


3  (x,  t)  u  (x,  t)  d¥- 


V 


f 


(t) 


)d¥  +  T ( n , x , t ) dS 


(t) 


sr  (t) 


(111) 


(112) 


and 


—  J  p(x,t)  { e  ( x , t ) 

Vf  (t) 


!  U(X,t) 


: d¥  =  j  :  (x,t)G(x,t).  u(x,t)d¥ 
vr  (t) 


+  T  (n  ,x  ,  t)  •  u  (x  ,t)  dS  +  j  q-n  dS 


(113) 


Sf  (t) 


sf  (t) 


where  e  and  G  are  the  internal  energy  and  body  force  per  unit  mass,  T  is 
the  stress  on  an  element  of  surface  dS  with  unit  normal  n,  and  q  is  the 
flux  of  energy  through  the  surface,  for  example,  heat  flux.  The  integra¬ 
tions  are  carried  out  over  Vf(t),  S^(t),  where  the  superscript  indicates 
that  the  control  volume  moves  with  the  fluid.  We  notice  that  all  the  terms 
contributing  to  the  balance  of  any  of  the  conserved  quantities  are  volume 
or  surface  integrals. 

In  the  case  of  an  inviscid  fluid 


T(n,x,t)  =  -  p(x,t)n. 


(114) 
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The  surface  integrals  I  TdS  and  jf-uds  reduce  then  to  i  pnds  and  /  pu-ndS 


f  r 

which  yield J  grad  pd¥  andj  div  (pu)dV,  respectively,  when  we  apply  the 


divergence  theorem. 

Recalling  Reynold's  transport  theorem 


d  /  -*■  3v 

^JX(x,t)d¥=j  d¥  +  j  x  u  •  ndS 


(115) 


V*(t) 


V*(t) 


* 

S  (t) 


where  V  (t)  is  a  control  volume  whose  surface  elements  dS  move  with  arbitrary 
velocity  u  .  Notice  that  the  two  integrals  on  the  RHS  are  over  space  and 
therefore  depend  only  on  the  instantaneous  position  of  the  control  volume. 

Consequently,  the  integration  can  be  carried  out  over  any  control  volume 

★ 

which  happens  to  coincide  with  ¥  at  this  instant,  whether  it  is  fixed 
or  moving  with  another  velocity.  Denoting  the  fluid  velocity  by  u  and 

-►q 

the  control  surface  velocity  by  u  wf;  get,  using  Eq.  (115) 

r 


d_ 

dt 


r 

3  \ 

X  d¥  = 

v—  dV  + 

J 

Xu' • ndS 


dt"  XdV  =  j  ft  dV  +  J  >;lj9'ndS 


q  f 

If  ¥  coincide  with  ¥  at  time  t,  we  get 


(116) 


CT  xd¥ 


d  ->  f  ->-a  ■> 

~  x^¥  +  (u  -  UM)-ndS 

dt  J  J 


_d_ 

dt 


(117) 


t 
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When  G  =  0  and  q  =  0,  Eqs.  (Ill)  through  (113)  become 


(118) 


pu  dV  + 


r  ,-*-f  -+g '  , 

)U  [  (U  -  U  )  T.jdS  = 


pndS ,- 


(119) 


d_ 

dt 


o  (e  + 


-vf ,  2 
u  1 

— + 


i  -*-f  i  2 

u  -*-f 

+  1  -  )  (U 


s 


9 


s9) 


:.dS 


(120) 


Using  the  divergence  theorem  we  can  get  the  differencial  form  of  the 
conservation  equations.  However,  it  is  far  more  convenient  to  use  the 
integral  form,  because  a  numerical  scheme  based  on  the  integral  form  is 
already  conservative,  since  the  fluxes  leaving  one  control  volume  have 
to  enter  an  adjacent  one,  and  discontinuities  can  be  propagated  in  principle 
without  any  smoothing,  since  one  can  always  integrate  a  profile  including 
a  discontinuity,  in  contrast  with  differentiation.  Consider  Fig.  10, 
representing  a  uniform  fixed  one-dimensional  grid  and  a  continuous  density 
profile  incorporating  one  discontinuity.  If  we  knew  the  mass  in  the  hatched 
cell  and  the  velocity  at  interfaces  A  and  B,  Eq.  (118)  will  give  us  the 
rate  of  change  of  mass  within  the  cell,  and  hence  the  mass  itself  after 
an  infinitesimal  time  6t.  But  we  have  to  get  the  density  at  A  and  3  and 
the  velocity  for  the  next  time  step.  We  must  have  recourse  then  to 
"averaging"  procedures  to  get  the  density  from  a  known  cell  mass  and  "inter¬ 
polation"  procedures  to  get  the  values  of  the  interfaces  from  the  cell 
average  values.  Through  these  two  procedures,  errors  are  introduced. 
Finally,  we  have  to  use  a  finite  grid  in  any  case. 
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Equations  (118)  to  (120)  can  be  written  in  a  reduced  form  as 


d  *  *  ->-f  -*■  q  -*•  *  * 

—  ptW  +  p  (u  -  uy) • ndS  =  -  T  dS  +  G  d¥ 
at 


*  *  ~bf 

where  p  is  a  generalized  density  (p  denotes  p,  pu  and  E 


=  P  (e  + 


rf  i2 

u  I 


in  Eqs .  (118)  to  (120) ,  respectively) ,  T  is  a  generalized  surface  stress 

(T  =  0,  pn,  pu  -n)  ,  while  G  denotes  a  generalized  body  force  (G  =  0  in 
Eqs.  (118)  to  (120)).  The  two  integrals  on  the  RHS  are  referred  to  as 
source  terms. 

A  naive  "finite-integral"  form  solution  can  be  written  as 


“mass  within 

mass  within 

net  outgoing 

control  volume 

control  volume 

_ 

mass  flux  through 

+  [source  terms 

at  t  +  <5t 

_at  t 

-control  surface  — 

As  will  be  explained  next,  the  above  formula  is  supplemented  with  diffusion 
flux  terms  (actually  diffusion  and  antidiffusion)  to  improve  its  accuracy. 


ACCURACY 

The  above  mathematical  analysis  was  carried  out  assuming  a  fixed 

uniform  grid  and  -rr  +  u  =  0,  where  u  5  constant.  We  notice  also  the 
dt  o  9x  o 

absence  of  any  source  term  (inhomogeneous  part  of  a  conservation  equation) . 
The  analytical  solution  was  found  out  to  be 


n+1 

P 


AP 


n 


“  i  Sc 

where  A  =  e  ,  then  was  expanded  to  get  a  numerical  scheme  that  matches 
it  up  to  a  prescribed  order  of  S.  In  this  context  the  numerical  scheme 
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is  an  approximate  solution  of  the  whole  PDE,  in  contrast  to  schemes  which 

3  9 

approximate  - — alone  by  a  finite  difference  and  - —  alone.  By  getting 

OX  ot 

a  solution  of  the  PDE  as  a  whole,  we  mix  the  time  and  space  derivatives 
for  a  higher  order  scheme.  To  see  that,  let  us  expand  o  (t  +  ot,x)  in 
a  Taylor  series: 


0  (t  +  <St,x)  =  p  (t)  + 


3o  i  <5t 

o  t  —  t  — 
3t'x  2! 


„2 

3  P  ] 

2  ' 

3t  x 


(121) 


From  the  PDE 


3p  _  _3p_ 

3t  Uo  3x 


(122a) 


A  .  ±t±)  ,  X,.a  la,  . 

2  3t  3t  3t  o  3x 
dt 


3  ,3p  ,  2  32p 

-  u  —  (- — )  -  u  - r 

o  3x  3t  o  2 
3x 


(122b) 


Substituting  Eqs.  (122)  into  (121)  we  obtain 

u  25t2  ,2 

o  (t  +  6t,x)  =P(t)  -uo5t^^-ii  + 

3x 


(123) 


showing  that  we  can  get  a  better  solution  in  the  time  domain  (of  higher 

3d 

order  in  6t)  by  adding  to  [o(t)  -  u  it  — ]  a  diffusion  term, 

o  3x 

2.2  go 

u  ot  .2  u^St) 

O  i  £  O 

- - - —  ,  a  purely  spatial  derivative.  Notice  that  - - -  is  equivalent 

3x" 

to  q-  ,  the  coefficient  of  dpn  in  the  schemes  discussed  earlier.  The  re- 

3  3 2 

maining  terms  appear  when  we  try  to  express  and  — ~  in  terms  of  finite 

dX 

3x 


differences  accurately. 

A  scheme  which  splits  the  time  and  space  domains,  on  the  other 
hand,  treats  Eq.  (122a)  as  an  ODE,  where  the  right-hand  side  is  assumed 
to  be  a  function  of  time.  A  second-order  Runge-Kutta  explicit  scheme 
can  be  written  as 
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p  (t  +  6t,x)  =  p  (t,x)  -  u  St  ~ !  3t  (124) 

O  aX  t  +  2  /  X 

where  —I  3t  is  obtained  by  first  getting  a  provisional  value  of 

oX  t  +  ~~  ,  X 

2 

the  density  at  t  +  ot/2  using  a  lower-order  scheme 

0  (t  +  4r>x)  =  p(t,x)  -  u  ~  irH  (125) 

2  O  2  a  X  t  ,X 

then  getting  3t  by  differentiating  p(t  +  ~  ,x)  spatially,  with 

9x  t  +  ,  x  - 

the  result 

—  I  dt  =  i£|  -  u  —  —I 

3x  't  +  —  ,x  3x't,x  o2  „2' 

2  ox  t ,  x 

Upon  substituting  in  Eq.  (124) ,  this  yields  Eq.  (123)  again.  One  can 
deduce  therefore  that  up  to  a  given  order,  schemes  which  mix  the  time 
and  space  domains  and  those  which  split  them  are  equivalent.  A  warning, 
however,  is  in  order  here:  A  concept  derived  for  a  split  time-space  scheme 
cannot  be  applied  directly  to  one  that  mixes  both  domains.  For  example, 
using  a  half  point  density  in  Eq.  (123),  i.e.,  the  scheme 


L  (t  ,  .  6t  3o  |  o 

P(t  +  T,X)  =  p(t,X)  -  Uq  —  -|t<x  +  ~g 


2  2 

u  Z6t  2 

O  dO 


(126a) 


3x  t ,  x 


p  (t  +  <5t,x)  =  p(t,x)  -  u  6t 


u  5t  ,2 
o  a  ; 


o  3x  t  + 


?x-  t+f  ,x 


(126b) 


will  cause  a  decrease  in  accuracy  instead  of  improving  it,  as  can  be  seen 
from  differentiating  Eq.  (126a)  with  respect  to  x  and  substituting  in 
Eq.  (126b) .  The  key  point  is  that  Eq.  (123)  is  a  solution  of  the  PDE 
as  a  whole. 


In  summary,  the  schemes  derived  in  earlier  sections  are  solutions 

f  q 

of  the  conservation  equations  if  u  =  constant,  u  =  0,  and  source  terms  =  0. 
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If  these  are  not  satisfied,  a  correction  that  preserves  the  order  of  the 
scheme  should  be  adopted.  Here  we  split  the  two  effects: 

(1)  ur  variable  and  source  terms  are  variable  j-  0 

(2)  u9  *  0 

and  treat  each  separately. 


GRID  MOTION 

According  to  the  above  splitting,  we  need  to  consider  a  case  where 
f  g 

u  =0  and  source  terms  vanish,  but  u  ^  0.  This  is  a  static  field,  where 
the  density  and  energy  are  constants.  Equations  (118)  and  (120)  reduce 
then  to 


This  exhibits  the  formula  for  an  accurate  scheme  when  the  grid  is  moving: 

the  rate  of  change  of  volume  equals  the  rate  of  sweeping  by  the  moving 

surface,  as  illustrated  in  Fig.  11.  Here  we  can  achieve  an  infinite- 

mean 

order  accuracy  in  ot  by  defining  a  mean  control  area  S  such  that 


smean 


t-ndS  =  swept  volume 


Let  us  consider  the  three  cases  of  1-D  geometry;  namely,  planar,  cylindrical, 
and  spherical  symmetries,  denoted  from  now  on  by  a  =  1,  2,  and  3,  respec¬ 
tively. 

In  the  planar  case,  the  area  is  independent  of  the  radius,  so 


that 
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=  1. 


where  L  and  R  denote  the  left  and  right  interfaces  of  the  control  cell 
respectively. 

In  cylindrical  1-D  geometry,  the  volume  swept  by  the  interface 


is 


r  ,  n+1  2  .  n  2. 

A¥B  =  Tt  (rB  }  "  (rB)  1  ' 

where  B  indicates  L  or  R.  Here  the  depth  of  the  cell  being  considered 


is  taken  equal  to  unity.  Since  5t  =  r^+1  -  r^,  the  average  area  is 

Q  B 


^n+i  _ 

B  "  ug5t 


r  .  n+1, 2  ,  n,  2. 

rB  ~  rB  1  ,  n+1  n, 

n+1 - n -  =1T(rB  +  V 


B 


-  r 


B 


One  can  define  then  average  radii 


n+i  _  1.  n  n+1, 
rB  =  2  ,rB  *  rB  >  ' 


„n, n+i, n+1  _  n, n+i, n+1 

since  A_  =  2+r 

B  B 


Finally,  in  spherical  geometry,  the  swept  volume  is 


4  n+1  3  ,  n,  3, 

=  -  tt[  (rQ  )  -  (rQ)  ]  , 


yielding 


.n+i  _  A¥  4  r ,  n,  2  .  n,  ,  n+1.  .  n+1. 2, 

'b  =  H+I - +  (rB,(rB  >  *  !l 


'B 


B 


whence 


n+i  ,1  n,2  ,  n, ,  n+1,  n+l,2,,i 

rB  =  3"  rB  +  (rB)(rB  5  +  (rB  )  ]}  ' 


since 


(128) 


(129) 


(130) 


(131) 


(132) 
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n ,n+i ,n+i  n,n+l,n+l  2 

AB  =  4  7t(rB  1  * 

Equations  (128)  ,  (129) ,  and  (131)  should  be  used  as  the  proper  interface 
areas  when  evaluating  the  fluxes  and  surface  forces.  To  complete  the 
formulation,  when  body  forces  are  present,  the  volume  used  should  be  that 
confined  between  the  average  interfaces.  It  can  be  arbitrarily  selected 
for  a  =  1,  and  is  defined  as 


v"+1  =rf(r"H)2  -  (r^*)2] 


for  x  =  2 ,  and 


»  §  .Ur"**)3  -  (r"+Vl 


(133) 


(134) 


for  a  =  3.  This  choice  will  ensure  a  proper  balance  between  surface  and 
body  forces. 


Variable  Velocity  Field  and  Source  Terms 

To  account  for  these  two  effects,  the  fluid  velocity  and  source 
terms  used  in  the  "finite-integral"  solution  should  be  evaluated  at  some 
intermediate  time  between  t“ ,  tn  ^  so  as  to  preserve  the  accuracy  of  the 
scheme.  Since  we  split  the  effects  of  grid  motion,  variable  velocity 
field  and  source  terms,  the  above  intermediate  values  should  be  derived 
from  an  ODE  solver  of  a  consistent  order  in  5t.  For  a  fourth-order 
(diffusion  and  phase  error)  accurate  scheme,  for  example,  we  need  a  second- 
order-accurate  explicit  ODE  solver.  In  other  words,  for  the  system  of 

-hf 

Eqs.  (118)  to  (120)  ,  we  advance  the  time  one-half  step  using  u  =  u 

n  .  n+i  .  -*f.n+i  _n+i  rT  ,  -+n+J  _  ,  -*-f^n+i,  n+} 

and  p  =  p  to  get  p  ,  (ou  )  ,  E  .We  define  u  =  (pu  )  /p 

and  pn+^  3  p(pn+^,  en+*)  where  p(p,e)  is  the  equation  cf  state  and 


i  til  "‘■n+i  ]  2 

en+i  =  (En+i/pn+i)  -  ^—-L 


-+£ 

Then  we  advance  the  system  a  whole  time  step  using  u  =  u  and 
p  =  pr‘+t  As  explained  earlier,  we  need  not  and  should  not  update  p, 
z u  and  E,  during  the  full  time  step,  since  the  scheme  is  already  a  solu¬ 
tion  of  the  whole  PDE.  For  a  sixth-order-accurate  scheme,  we  need  a 
fourth-order  ODE  solver  and  so  on. 


Example  of  an  Algorithm 


Let  us  implement  the  scheme 


n  e  ,  n  n  .  ,  £  ,  n  nn 

:3  -  I(pj+1  -  Dj-1)  +  4  (pj+r  2pj  +Pj-1); 


TD  T  1  e  ,  ,  n  „n  n, 

+  (6  +  -  2Dj  + 


n+1  TD  1  £  ,  T  T  T  , 

Sj  ‘  °j  -  <6  -  6  1<Cj+l  -  20i  *  11 

a  stable,  fourth-order  phase  error,  sixth-order  diffusion  error  scheme, 

where  p denotes  either  of  p,  ou  ,  or  E. 

If  we  have  N  cells  whose  interfaces  are  at  radii  (rP  ,  r1?^,..., 

1  /  -  3/ 

a  ,  .  n  .  ,  n+1  n+1  ,  ^  .n+1  ,  . 

r  .  at  time  t  ,  moving  to  (r,  r.  )  at  t  ,  let  us  denote 

N+l/2  1/2  3/2 

the  cell  centers  by  the  subscripts  j  =  1,  2,...,N,  located  at 


n  ,n+l 


1  ,  n ,n+l  n ,n+l,  .  ,  _ 

2<rj-l/2  +  rj+l/2  £°r  “  *  ‘'2 


,lr/  n,n+l,2  ,  n,n+l.  ,  n,n+l.  ,  n,n+l,2,,J  , 

J  rj-l/2  *  'I]-l/2Ucj*l/2)  *  (t  3*1/2'  ]}  f°r  '»  '  3- 


The  volume  of  the  j  cell  per  unit  angle  is  given  by 
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n+1  T  n  n  ,  n  n+i  ..  n+i  rM_.  n  n+i  _  n+i. 

V .  o  .  =  V .  D  .  -  ot  (p  .  ,  A .  .  6U  .  f )  +  5 1  (o  .  .  A  .  ,  OU .  : ) 

3  3  33  3  +  i  3  +  i  3  +  i  3“i  3~i  3"i 


.  n+i  .n  ,  n  n,  ,  n+i  ,  n  ,  n  n  n+i 

+  A  V  .  (p  ,  -  p  . )  -  X  .  ,  V .  , (p .  -  p .  )  +  source  . 

3  +  i  3  +  i  3  +  i  3  3“i  3-i  3  3"1  3- 


where 


n  1  .  n  n 

•Vi  'i<aj  *  Vi 


for  j  =  1,...,  N-l,  while 


n  1  .  n  n, 

:i  =  2(?L  +  °l)  ' 


"N-i 


1  .  n  ,  n. 

—  (p  +  p  )  , 

2V  N  r'  ' 


where  L  and  R  denote  left  and  right  guard  cells,  respectively.  The 
difference  5U_.  +  ^  between  the  fluid  and  grid  velocities  is  given  by 


...n+i  ,,n+i  -a.  c a.  tTn+i  tA  ,  n+1  n  , 

OU.  ,  ot  ~  U.  ,  ot  -  U.  .fit  =  U.  ,  5t  -  (r.  ,  -  r.  .), 
3+i  3+i  3+i  3+i  3+i  3+i 


while  the  diffusion  coefficient  is 


Xn+i  =  i(£n+ * ) 2 
3  +  i  4  G  j  +  i  ' 


where 


.n+i 

"j  +  i 


.  n+i  n+i  . 
5U  .  , i  A.  ,  Ot 


&  —  (JL  + 


vn  vn  , 

3  3  +  1 


The  velocity  at  the  interfaces  satisfies 


«■:!  ■  K‘  *  -#> 


for  j  =  1,...,  N-l,  while 


(141) 


(142) 


(143) 


(144) 


(145) 


(146) 


«  - 
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n+i  n+i  n+i  n+i 

V  =  UL  '  Vi  =  UR  ■ 


The  volumes  V.  .  are  defined  as 
3  +  i 


vn 


3  +  i 


■  *  ^x> 


(147) 


for  j  =  1,...,  N-l ,  and 


Vn  =  ¥n ,  V  =  ¥“  . 
i  l'  N+i  N 


Equation  (135b)  then  adds  the  main  diffusion,  giving 


„n+l  TD  n+1  T 

V .  p  .  =  V  .  p  . 

3  3  3  3 


n+i  ,,n 


n+i  ,,n 


+  v .  :  ¥.  , (p .  ,  -  p .)  -  v .  :  v.  (p . 
3  +  i  3+i  3  +  -’  3  3-i  3~  i  3 


h-i> 


(148) 


where 


n+i 

Vi 


(rn+i)2 

I  +  -  j±f_ 

6  12 


(149) 


Finally,  the  antidif fusive  fluxes  are  evaluated  according  to 


_  n+i  n+1.  T  T, 

r  3+i  Uj+i  j+i  ° j  +  i  K j  ' 


(150) 


where 


n+i 

Vi 


1_ 

6 


,.n+i  2 
+j + i 


(151) 


and  then  corrected  using  one  of  the  flux  limiters  Eq.  (64)  or  Eas.  (65) -(75) 
Let  us  select  Eq.  (64)  on  account  of  its  simplicity.  The  corrected  fluxes 
are  given  by 

c  n+1  TD  TD 

F.  =  sign  (F .  , ) • max  {0,  min  [sign  (F.  ,)■¥'•  (p  -  p  ), 

3  +  i  3  +  i  3+i  3  +  i  3++  3+1 


I F  .  .  !  ,  sign  (F  ,  ,  ) • v .  , 

1  3+J  3+i  3-i  3 


>¥n+l  (oTD_  oTD  ,, 


5-1" 


(152) 


whence 


wr 
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n+i  _  ^TD 


n+1  (F]  +  i  ‘ 


(152) 


As  for  the  source  terms,  they  are  summations  over  the  surface  cr  the  volume 
of  the  cell.  Lot  us  consider  first  l-j pn  dS),  which  yields  [-  grad  p] .  In 

5 

cartesian  coordinates,  following  the  diagram  of  Fig.  II, 


n+i  n+J  n+i  n+i  n+i 

source  =  d  ,  A .  .  -  p  .  ,  A .  , 

3  3-1  3“i  3  +  i  3  +  i 


(154) 


where 


n+i  1 ,  n+i  n+i 

pvt  ■  l<pj  '  EVil 


(155) 


,  .  „  .  ...  n+i  n+i  n+i  n+i 

for  3  =  1,...,  N-l,  while  p^  =  pr  and  p^  =  pR 

In  cylindrical  geometry,  following  Fig.  13,  we  have 


n+i  n+i  n+i  n+i  .n+i  n+i  ,  n+i 

source  .  =  p.  .A.  ,  -  p.  .A,  .  +  p.  (r.  . 

3  3-i  3“i  3+i  3+i  3  3+i 


and  since 


n+i  n+i 

r .  ,  -  r  .  . 
3+i  3-i 


n+i  2  n+i  : 

rn:|  * 

3+i  3~i 


V1 

3 


+i 


n+i' 


where  from  Eq.  (136) 


n+i  1,  n+i  n+i, 

r  •  =  ~(r  *  +  r .  .) 

3  2  ]+i  3-i 


(156) 


for  a  =  2,  we  can  rewrite  the  expression  for  source  as 


n+i  n+i  n+i  n+i  n+i  ^ 


.n+i 


source  .  =  p.  ;  A'."!  -  p'."!  A"  f  +  — 2 — -  ¥n+‘ , 

3  3-i  3-i  3+i  3+i  rn+i  3 

j 


(157) 
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are  obtained  by  first  advancing  the 


Finally,  un+^  and  source  n+^ 

3  3 

whole  system  of  Eqs.  (136)-(160)  a  half  time  step  using  U^,  source  then 
a  whole  time  step  using 


0n+i|  =  ^i+l, 

j  ‘t-*-t+5t  j  lt->-t+6t/2 


(161) 


and 


source 


n+i  | 

j  1  t-*-t+6t 


=  source 


n+1 1 


'  t-»-t+6t/2  ‘ 


(162) 


I 


XI.  TWO-DIMENSIONAL  TRANSPORT 


Now  let  us  consider  the  two-dimensional  equivalent  of  Eq.  (2) , 


3p  9p  3o  n 

3t  lSx  2W  °' 


(2C1) 


whose  analytic  solution  is 


p(x,y,t)  =  p(x  -  uxt,y  -  u?t,Q) , 


(202) 


a  wave  propagating  with  velocity  u  =  (u  ,  u,) .  Assuming  an  initial  density 

1  C. 

c(x,y,0)  =  F(x,y),  we  Fourier  analyze  F(x,y)  in  space  on  a  rectangle  x  L0 
with  periodic  boundary  conditions: 


F(r) 


.n. 

k=-  ® 


-ik  ■ 


(203a) 


where  r  =  (x,y),  and  k  =  (k  ,k,,)  is  assumed  to  be  normalized,  i.e.,  k 
kl  k2 

denotes  2 tt  ( — ,  — )  .  Notice  that  the  summation  of  (203a)  is  actually  a 
L1  L2 
double  summation. 


- - ,00  r - ®j 

F(X'y>=L  L  Vk: 


i  (k  x  +  k^y) 
el  t- 


(203b) 


Lr  o—co  V  ——co 
1  2 


To  gain  insight,  let  us  consider  only  one  wave  component  of  Eq. 


(203)  , 


F(r)  =  sin  k-r 


(204a) 


or 


kjX  k^y 

F  ( x  , y )  =  six  2 7t  (- -  +  - - ). 

L1  L2 


(204b) 


Figure  21  shows  the  resulting  waves  for  different  values  of  (L^L^),  (k^,k.J. 
From  Eq.  (204b),  F(x,y)  is  constant  along  lines  of  constant 
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klX  k2y 

(- —  +  - - ).  For  example,  the  nodes  of  the  wave  coincide  with  the  lines 

L1  L2 


V  ^2*  1  .  3  , 

J  -N  /  *-  /  ^  /  Z  t  m  •  m 

L1  L2  - 


(205) 


which  are  normal  to  the  wave  vector  k. 


To  find  the  wave  length  for  a  given  (k^/k^),  we  first  go  back  to 

the  one-dimensional  case.  For  a  system  of  length  L  and  periodic  boundary 

kx  kv 

conditions,  the  harmonics  sin  2^-—  and  cos  2tt-t are  admitted,  where  k  =  0, 

La  Li 

1,2,...,=°.  With  each  of  these  is  associated  a  wave  length  X  defined  as  the 

kx  27ikx 

distance  between  two  successive  "even"  nodes.  Since  sin  2*r— —  =  0  at  -  =  0, 

L  Li 


m ,  2 m  ,  3 m .  \  is  obtained  from 


2irkX 


2tt,  yielding 


'4- 


(206) 


We  get  therefore  wave  lengths  *>,  L,  L,  Lf  . . .  ,  where  the  longest  finite 
wave  length  equals  L,  the  system  length.  In  two-dimensional,  the  wave 
length  for  a  given  k  is  defined  analogously  as  the  distance  between  two 

—V 

points  on  successive  "even"  node  lines,  projected  on  the  direction  of  k. 
From  Eq.  (204b) , 


F(x,0)  =  sin  2tt 


k^x 


which,  as  explained  above,  yields  X  =  —  where  A  is  the  wave  length 

1  k!  >S 

along  the  x-direction,  which  when  projected  on  k  =  2tt( — ,  -=-)  yields  \: 

L1  L2 

L1  kl  k2 

(r-,0).(-^,  -^) 

kl  L1  L2 


X  = 


2  tt 

Ikl 


(207a) 


*  S. 
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where 


k 


2Tk,  2itk 

19  9  9 

< T1'' 

1  2 


or 


A  = 


.  1,2  12 

(t~)  +  (t-> 

x  Av 


(207b) 


^1  ^2 

N  w  we  find  all  the  wavelengths  along  a  given  direction  - —  :  —  = 

L1  L2 

constant  c.  Noticing  that  ,  k^  for  periodic  boundary  conditions  can  take 
only  integer  values,  the  waves  along  a  given  direction  correspond  to 


(n) 


.  (1)  .  (n) 

nk1  , 


nk 


(1) 


(1)  ,  (1) 


^  (n  =  1 ,2,...,co)  where  k^  ,  k?  are  the 

smallest  integers  that  satisfy 


Prom  Eq.  (207a) 


— 

fnk  ,  nk 

— >“  +  T 


(208) 


wnere 


k(1)  k(1) 

1  2  2 
(-7—)  +  (~ -)' 


A  X_i 

1’  2  '  3  " 


A, 


(209) 


1  2 

Along  a  given  direction  we  have  wave  lengths 
Consider,  for  example.  Fig.  21(b),  where  L1  *  2 ,  =  1.  Along  direction 


(1/2,  1),  k^  =  =  1,  whence  A1  =  1  /[  1  2  +  =  2/\Js.  The 


maximum 


system  length  along  this  direction 


,  .  ^[J7 2  7  \R  X1  2,\£*  4 

being  +  i  =  =  _ 


LJ,1  ^*/2  5 

showing  that  because  of  the  periodic  boundary  condition  independently  in 


each  direction  the  longest  wave  length  is  only  80  percent  of  the  maximum 
system  length  in  the  direction  (-j,l),  in  contrast  to  one-dimensional 
cases  where  A^  =  L.  For  the  case  of  Fig.  21(a),  A^  is  50  percent  of  the 
system  length. 

From  Eq.  (202) , 


0  (} 


.  .  \  x  i[k,  (x  -  u.t)  +  k  (y  -  u.t)  ] 

x’y't>  ■  L  L  \.k  ’  1  1  2  2 


k  =— co 

l  k2»- « 


-E  E 

ki  k2 


Pk  k(t)  ei(klX  +  k2y) 
1'  2 


(210a) 


or 


-E 


ik-  r 


k=- 


where 


(210b) 


'  ...  _  "  -i  (k  u  +  k  u  )t 

pv  v  e  11  22 

/*2  ^ 


(211a) 


or 


P^(t)  =  e 
k  k 


-  ik  •  ut 


(211b) 


Thus  each  harmonic  independently  advances  uniformly  in  phase  without 
changing  its  magnitude,  as  shown  in  Fig.  (22). 
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We  notice  that  the  different  harmonics  advance  in  the  direction 
u,  which  is  generally  different  from  that  of  k,  as  illustrated  in  Fig. 

(23).  They  keep  their  front  normal  to  k  and  therefore  the  projection 
cf  u  on  k  is  the  speed  cf  advance.  This  adds  extra  requirements  that 
were  not  invoked  in  the  one-dimensional  case,  namely: 

1.  The  scheme  should  kern  the  wave  front  a  straight  line; 
otherwise  distortion  of  profiles  occurs. 

—f" 

2.  It  should  also  keep  the  wave  front  normal  to  k;  otherwise 
"scattering"  occurs,  namely  waves  with  different  !k!  but 
the  same  direction  (k^:k^)  will  come  out  in  different 
directions,  causing  scattering  of  the  transported  profile. 

As  will  be  proved  later,  the  speed  of  propagation  V  of  a  numerical 
scheme  differs  from  u  not  only  in  magnitude  as  in  one  dimension  but  also 
in  direction,  providing  one  more  source  of  error.  If  the  above  two  require- 

->  -V  -> 

ments  are  satisfied,  however,  only  k* (V  -  u)  contributes  to  the  phase 
error. 


Now  supoose  is  known  at  all  times  only  on  a  set  of  (N.  +■  1}  •  (Nn  t  1 

L,  L, 


discrete  grid  points  with  separation  6x  =  — ,  (y  =  namely,  x.  =  i'x 

N  ^  ^  1 


(i  =0,1,...,  N^-  1),  y  =  j?y  (j  =  0, 1 , . .  .  N?  -  1),  the  origin  being  a  member 
of  the  set.  According  to  periodicity  assumption  c 


o,  ■ 


N1H2 

hence  we  can  have  only  — - —  +  1  different  harmonics.  Let 


N  , j ,  i,o= 

1  l,^.; 


N, 


N. 


f  ( 


.  2  r-1  2  5  i  (k.  x  + 

lX'y>  ’L  -»aL  1 


k2y) 


(212) 


kl  = 


k2  =  7 
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Since 


i(k.x  +  k.,y ) 
e  1 


N  21k,  k 

2^!  Nl  -  ei(q(T)i5x  3'7)  =  ei<-+2,-  joy) 

ki  =  T 


2  .  . 


2irk, 


i[iir  +  2 tt— - —  joy)  -  2iri]  _  gi  (-itr 


j5y) 


-N  27Tk 

=  ei[r(~)i5x  +  ~  j6y]  =  ei(klx  +  k2y)  '  "1 

1  2  ki  =  T 


-N, 


for  all  i,  and  similarly 


ei(k1x  +  k2y)  |  -N2  _  gilk^  +  k2y)  |  N2 

k2  =  I"  k2  =  T 

for  all  j,  Eq.  (212)  can  be  rewritten  as 


N, 


f  (x,: 


y)  *E  2  £ 

-N 


*  i(v^2y) 

kl'k2 


-N, 


ki=_T~  +  1  V  ~T  +  1 


showing  that  k  space  structure  contains  only  x  N2  independent  points 

(see  Fig.  24).  The  amplitudes  6  can  be  obtained  from 

kl'  2 


•E-E2  E 


-N  -N 

kr  ~ + 1  V  —  + 1 


2irk  i£x  2’tk  j5y 
1  i  (— -  +  — r - ) 

W  1 


(214) 


for  i  =  0,1,2,...,  N^-l,  and  j  =  0,1,2,...,  N0-l. 


In  terms  of  sines  and  cosines,  Eq.  (213)  can  be  written  as 


f =  Ao,o  + 


_1  ,  _2 

E 2  V 

V1  V1 


2^TK1X  2TTKn  y 

cos  (  .  —  +  r  > 
ko  h  l2 


2irk1  x  2irk0y  2ffkix  27Tk2y 

*  \  k  sin  (Ti  *  -r-)  *  Ck  «»(— - — > 

kl,K2  1  2  v  12  /I  z 


N  / 

2iTk,x  2rtk7y  >  T  "  1 1  ,z1TklX 

* v*2  ■“  <- h5-  -  tt*  rL  V°  ~ 

12  \  /  V1  \  ; 


2ifk  x  ^x  r^1!- 

+  B  (sin— - )  +  AN  "oslE  )  +  /L/ 

kl'°  Li  /  _i0  1  v  -i 

/  ~'°  k2  X 


2irk  y  ™2y 

,,k  ain  'TJ-’  *  Vn2cos(  l2  h\  n, 

2  2  '  2 


A  ,  cos  (- 
o,k? 


where 


cos  (- 


kl'k2 


\  ,k_  "  iBk  ,k, 

_ i 2 - i — -  for  k,  >  0,  k  >  0 

2  12 


Ak.  ,k.7  Bk,,k 

_ l — - = — ~  for  k.  <  C,  k  <  0 

2  12 


C,  ,  -  iD,  v 

_J: i - = — ~  for  k.  >  0,  k  <  0 

2  12 


Ci<^  ,k„  +  lDk1<k-, 
- 


for  k  <  0,  k2  >  0 


■  ,  -  -  •  •  r  ■ 


•Vb  * 
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Again,  we  have  in  Eq.  (215)  x  N7  coefficients:  A^  (k^ 


(k,  =  0,1,...,-^ 


and  k0  =  0,1, 


12 

N„  N  N, 

B  (k  =  0,1,...,— - 1  and  k  =  0,1,...,^-  1), 

Ki'K2  “  2 


N  N 

1  2 

Ck  k  '  °k  k  ^kl  =  ' ' 2 '  *  *  ‘  ,~2~  ~  1  and  k2  =  1  >2 '  •  ■  •  ~  that  can  be 

determined  from  the  system  of  equations  f(x. ,x.)  =  p°  .. 

i/l 

Going  back  to  Fig.  24  let  us  count  the  different  harmonics.  The 
harmonics  are  considered  equal  if  they  have  the  same  magnitude 

2Trkl  2  27rk2  2  kl  k2 

— )  +  (— . )  and  are  aligned,  i.e.,  - —  :  —  =  constant.  The 

L1  L2  L1  L2 

number  of  the  harmonics  is  almost  half  the  space  of  Fig.  24  since  (k. ,k  ) 

±2 

is  equivalent  to  (-k^-k^)  and  (k^-k^  is  equivalent  to  (-k^k^).  For 

example:  a  and  b  in  Fig.  24  are  equivalent.  Figure  25  shows  the 

independent  harmonics  selected  to  match  the  choice  in  Eq.  (215)  .  The 

N1N2 

number  of  the  harmonics  is  therefore,  — - —  +  1. 

If  we  count  the  maximum  number  of  wave  lengths,  we  get  an  even 

2  a 


smaller  number,  since  according  to  Eq.  (207a) ,  X  = 


Two  harmonics 


Ik 


such  as  a  and  b  in  Fig.  25  will  give  the  same  value  for  k 


The  maximum 


N. 


N„ 


number  of  wave  lengths  is  therefore  (—  +  l)My  +  I),  corresponding 

to  the  positive  quadrant  of  Fig.  25.  This  is  an  upper  limit.  This  is 

because  the  number  of  wave  lengths  can  be  less  if  the  ratio  cx/6y  is  a 

rational  number.  As  explained  above,  decomposition  in  two  directions 

puts  a  limit  on  the  longest  finite  wave  length  X  in  a  given  direction. 

Discretization,  on  the  other  hand,  puts  a  limit  on  the  shortest  wave  length 

in  a  given  direction  since  it  reduces  n  in  Eq.  (208) .  The  largest  value 

N1  N2  (1)  p) 

occurs  for  k^  =  —  ,  k^  =  — .  If  k^  ,  k^~  are  the  smallest  integers 


for  a  given  direction,  the  shortest  wave  length  along  this  direction 
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corresponds  to 


.  N  /2  N./2 

mm  1  2 

n  "  integer  (1)  '  (2)  ‘ 

1  K2 


Assuming  p(x, y,0)  =  f(x,y),  i.e.,  asstaming  the  density  in  between 

the  grid  points  values  p°  .  to  be  f(x,y),  Eq.  (202)  predicts  the  density 

1  >  3 

at  time  t  as 


■:V,t)  *  £ 


5  i[k. (x-u  t)  +  k  (y-u  t) ] 

Q ,  ,  G  i-  j.  Z  Z 

Vk2 


-N  -N 

+  1  k2=~  +1 


=  E  Ev 


(t)ei(kiX  +  k2y) 


kl  k2 


,  .  ...  *  -i(k,u.t  +  k-u  t) .  Since  we  are  oniv  concerned 

where  o  (t)  =  p  e  1  1  2  2 

K1'K2  K1'K2 

with  o(x^,y_.,t)  ,  let  x  =  x  ^  =  i6x,  y  =  y^  =  j6y.  We  then  get 


E  z 

“=5“  A  i(k-iox  -r  k-jcy)  .  If  the  time 

2  pk(t)e  1  2 


-N  -N, 

kl=—  +  1  k2~  +  1 


is  also  discretized,  let  tn  =  not,  oa  .  =  p(x. ,y.,tn)  and  on  ,  =  p.  .  (tn) , 

ig  i  3  x,  ,k  k  ,k? 


Then  p 


n  =  p  m  i  (k  i6  x  +  k,; 

i,j  4-/  L-,  klfk. 


kl  k2 


where 


kl'k2  "kl'k 


i  (k  u  +  k  u  ) not 
e  II  2  2 


We  define  A(k^,k2)  as 


wmwwwm 


■Ve  * 


A(k  ,k  )  =  - 

12 


kl'k2 


Equation  (218)  expresses  the  analytic  solution  as 


A(krk2)  =  e"i(klUl5t  +  k2u26t)*  (23 

u^cSt  u^St 

If  we  denote  — - — -  by  e  ,  — ~ —  by  e  ,  k.6x  by  3  ,  and  k_<5y  bv  3  /  Eq.  (219) 
6x  x6y  yl  x  2  y 

reduces  to 

_  ,  -i  (e  6  +  e  3  )  (21 

A(kl'k2  =  6  X  X  y  Y 

Mow  let  us  analyze  a  fully  two-dimensional  scheme,  a  direct  extension  of 
the  one-dimensional  scheme 


T  n  e  ,  n  n  . 

■  h  -  i«y i  -  vi 


TO  T  ,  n  _  n  n 

°j  =  Pj  +  p(pj+l  ‘  2pj  +  ’j-l* 


n+1  TD  ,  T  T  T 

Pj  =  p.  -  W(P.+1  -  20.  +  0j_1) ' 


name ly , 


n  x  ,n  n 

Pi,j  2  Pi+l,j  °i-l,j 


n  Y  ,  n  n  . 

Pi,j  ~  2  Pi , j  +  1  -  °i,j-l)! 


n  x  ,  n  n 

o  .  .  — —  (o  .  .  -  o  .  . ) 

b]  2  i+1,d  i-l,D 


-Z  rnn 


IT  (P  ■  .  ,  i  ~  P  ,  4  ,  )  S 

2  1,3+1  1,J-1 


TD  T  ,  n  n  n.  ,  n  _  n 

P.  .  =  P.  .  +  v  (p.  ,  -  2o .  .  +  p.  .  .)  +  v,  (o .  .  .  -  2p  .  . 

1,3  i-3  x  1+1,3  1,3  1-1,3  y  1,3+1  1,3 


li  V 

+  Pi , j-1  ' 


TD 

i-3 
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n+1 

o .  . 

1-3 


,  Tx  „  Tx  Tx 

-  (p .  ,  .  -  2c  .  .  +  3  .  , 

x  1,3  1-3 


.) 


u 


y 


(p 


Ty 

i- j  + 


i 


^  Ty  Ty 

"“'i- j  ~i,  j-1)  . 


(222e ) 


Again  as  in  one-dimension,  after  r ourier-analyzing  the  initial 

density  profile,  i.e. ,  after  we  have  gotten  the  [3,  .  ],  the  problem  is 

,  l  ,K2 

i(}c  X  +  v) 

reduced  to  propagation  of  the  complex  harmonics  el  2'  .  Since  for 

the  linear  problem  Eg.  (201),  each  harmonic  propagates  independently,  we  can 
get  Afk^k^)  by  assuming  only  one  harmonic: 

n  _  o  i(k  i<5x  +  k  joy)  o  i  (i3  +  j3  ),  (223) 

o  =  e  1  2  =  o  e  x  v 

i- j 

then  using 

.n_rl 

K  • 

A  (k  ,k  )  =  .  (224) 

12  n 

0  .  . 
i-3 

Substituting  Eg.  (223)  into  (222a)  we  get 


Tx  o  i(k.i5x  +  k.jcy)  "x  ,  o  i[k.  (i+l)6x  +  k.j<Sy) 
o .  ,  =  p  e  1  2  -  — —  \  o  e  1  2 

i-3  2 


o  i[k.  (i-l)Sx  t  knj<5y)} 

”  p  0  i. 


(225  > 


r*ence 


Tx 

1 LI 

n 

) .  . 

1-3 


£  is  -iS 

x  v  X 

1  -  —  (e  '-e  ‘  )  =  1  -  i  e  sin  6 

2  xx 


(226a) 


Similarly  Eg.  (222b)  gives 


o^. 

i»3  ,  . 

- —  —  1  —  l  e  sm  3  . 


„  n 

"i- j 

Denoting  i  sin  8^  by  t^  and  i  sin  S^byt^,  Eg.  (222c)  gives 


(226b) 


—  ^  —  1  —  £  t  —  £  t 


n 

o .  . 

1-3 


xx  xx 


(226c) 
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Substituting  Eq.  (223)  into  (222d) 


TD  T  .  .  .  „  .  _  ■  „ 

P  .  .  P .  .  1 2  -id  ip  - iS 

— ibl  =  — i-tl  +  v  (e  x  -  2  +  e  x)  +  v  (e  y  -  2  +  e  y) 
n  n  x  y 

p .  .  p .  . 

1/3  1,3 


T 

P .  • 

=  -  ^  +  v  d  +  v  d  , 

n  xx  y  y 

P .  . 

1.3 


(226d) 


where  d^  =  2  (cos  6^  -  1)  and  d^_  =  2  (cos  -  1)  •  Finally,  Eq.  (222e) 

yields  with  Eqs.  (226c)  and  (226d) 

n+1 

P  ■  ■ 

A  (8  ,8  )  =  lj  =  (1  -  e  t  -  tt)  +  vd  +  vd  -  u  d  (1  -  -:t) 
xy  n  xx  yy  xxyy  xx  xx 

Pi.  j 


-  u  d  (1  -  £  t  ) 
y  y  y  y 


(226e) 


From  which 


S^+l.  =  A(8x,8y)p°el(l6x  +  jV 

i(iB  +  j21  ) 

We  notice  that  the  coefficient  of  e  x  vy  is  independent 

i  and  j,  i.e.,  independent  of  x,y.  Consequently,  pn  .  have  the  same 

1,3 


of 


wave  front  inclination  and  shape  as  f,  i.e.,  along  the  lines  of  k.x  +  k  v 

1  f  J  i  z 

=  constant,  pn+1  =  constant.  Finally,  this  is  a  nine-point  explicit  scheme, 
1.3 

as  illustrated  in  Fig.  26,  which  shows  the  points  involved  in  determining 

n+1 
P  . 

1.3 


- f  r  -1 


r 


i 


XII.  AMPLITUDE  AND  PHASE  ANALYSIS 

We  write  A  as  A  =  ;A|ell'/  where  J A |  is  the  amplitude  and  9  is 

the  phase  angle.  To  classify  the  order  of  the  scheme  we  need  to  expand 

| A (  and  6  in  a  power  series  in  3^  and  3  . 

PHASE  ERRORS 


In  the  two-dimensional  case  we  have  A  =  A (t  ,t  ,d  ,d  )  where  t  ,d 

x  y  x  y  xx 

are  functions  of  3  ,  while  t  ,d  are  functions  of  3  •  Since  loa  A  =  log  !a! 

x  y  y  y  •  ' 


+  id), 

3  =  In  (log  A)  . 


(229) 


Expanding  3  in  a  power  series  of  3  ,  3  near  3,3  =  0,  we  get 

x  y  x  y 

S2  3 2 

3  =  5  +  (3X3  +  3y3  )  +  ( 3XX  •—  +3Xy3  S  +  3  yy  ^  ) 

o  ox  oy  o2  oxy  o2 

2 3  ,2  2  3 

,oxxx  I*  ■  ,*xy  xpy  .xyy  x  y  ayyy  ,  , 

^  \  “  _  r  ~  ^  T  „  )  *r  ...  , 

o  6  o  2  o  2  06 

where  3X  =  (33/33  )  at  3  =0,  while  3y  =  (33/33  )  at  3  =0,  etc. 

o  xx  o  y  y 


(230) 


(230) 


We  therefore  need  the  derivatives  of  log  A.  Noticing  that 


'*y= 

0)  = 

1  we 

get 

(log 

„  X 

x 

A)  = 

A 

(231a) 

0 

o 

(log 

A)  y  = 

Ay; 

(231b) 

o 

o 

(log 

„  XX 

XX 

,  x.  2 

A)  = 

A  - 

(A  )  ; 

(232a) 

o 

o 

o 

(log  A)Xy= 

o 

XV 

A 

o 

AXAy; 
o  o 

(232b) 

A  =  t  t  A  1  +  t  d  A  +  d  t  A  +  d  d  A 

o  xo  yo  o  xo  yo  o  xo  yo  o  xo  yo  o 


III  tv  III  dy  I  II  tytV 

a  Jt  4.  -3  4-  4-  £  * 


Ayyy  =  t  A  y  +  d  A  y  +  3t  t  A 
o  yo  o  vo  o  yo  yo  o 


I  II  II  I  Cydy  '  "  dydy 

+  3 (t  d  +  t  d  )A  y  y  +  3d  d  A  y 


yo  yo  yo  yo  o 


yo  yo  o 


Finally , 


vvvv  'V  fcv  IV  dv  <  *"  11  V  tt  llt  , 

A  =tAX  +  dAX+(4tt  +  3t  ""  )  A  X  X  +  (4t  d 
o  xo  o  xo  o  xo  xo  xo  o  xo  x< 

II  II  I  III  t  d  I  ,11  111  dyd,, 

+  6t  d  +  4t  d  )A  +  (4d  d  +  3d  )A  ; 
xo  xo  xo  xo  o  xo  xo  xo  o 


,,,  ,  t  t  ,  t  d  d  t  d  d 

Axxxy  =  t-  (t-  A  x  y  +  d'  A  x  y,  +  d'"  {t  A  X  y  +  d  A  *  y 


XO  VO  o 


xo  yo  o 


xxvv  "  "  Vv  "  "  dxdv  "  "  txdv  "  " 

y  _  ^  +  d  d  A  *  +  t  d  a  ^  +  d  t  A 

o  xo  yo  o  xo  yo  o  xo  yo  o  xo  yo  o 

...  ,  t  t  ,  d  t  ...  .  t  d  ,  d  d 

A*YYy  «  t  (t  Axy+d  AXy)  +  d  (t  A  X  y  +  d  A  X  y ) 
o  yo  xo  o  xo  c  yo  xo  o  xo  o 


yyyy  1 V  fcy  •  V  dy  ■  "*  "2  W  „ 

Ay//y  =  t  A  y  +  d  A  y  +  (4t  t  +  3t  }A  y  y  (4t  d 

yo  o  yo  o  yo  yo  yo  o  —  yo  yo 


II  II  I  III  t-yy  I  HI  1*2  Qydy 

+  6t  d  +  4t  d  )A  y  y  +  (4d  d  +  3d  )A  y  y; 


yo  yo  yo  yo  o 


yo  yo  yo  o 


Going  back  to  the  definitions  of  t  ,t  ,d  ,  and  d  , 

x  y  x  y 


t  =  t  =  0; 
xo  yo 


t  =  t  =  i; 
xo  yo 


t  =  t  =  0; 
xo  yo 


lit  III 

t  =  t  =  -i; 
xo  yo 


IV  •  v 

t  =  t  =  0; 
xo  yo 


J.  A-  .  -A 


^  V;  * 
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d 

= 

d 

=  0; 

XO 

yo 

d 

1 

d  = 

=  0; 

XO 

yo 

If 

■  i 

d 

= 

d 

=  ~2; 

XO 

yo 

■  1 1 

■i  i 

d 

d 

=  0; 

XO 

yo 

'  V 

IV 

d 

= 

d 

=  2. 

XO 

yo 

(240b) 


t  t  t  t 

xx  y  y 

Substituting  into  Eqs.  (235) -(239)  and  assuming  A  =  A  =  0,  we  get 


A_  =  1; 


(241) 


o 

t 

t 

AX  = 

iA  X, 

v  y 

A'  =  lA 

(242) 

o 

o 

O  O 

d 

t  t  d 

AXX  = 

-2A  X, 

A^  =  -  A  X  y,  Ayy  =  -  2A  y; 

(243) 

o 

o 

o  o  o  o 

XXX 

A 

t  t  d 

=  -  i(A  +  6A  ^  ; 

(244a) 

o 

o  o 

t  d 

Axxy 

=  -  2i 

A  y  X; 

(244b) 

o 

o 

t  d 

AXYy 

=  -  2i 

A  X  y ; 

(244c) 

o 

o 

t  t  d 

Ayyy 

o 

=  -  i(A  y  +  6A  y  y); 
o  o 

(244d) 

and 


d  d  d 

xxxx  x  x  x. 

A  =  2 (A  +  6A  ) ; 
o  o  o 


(245a) 


t  t 

.XXXV  .  X  V 

A  *  =  A  J  ; 
o  o 


(245b) 


d  d 

AXxyy  =  4A  x  y 
o  o 


AXyyy  =  A 


t  t 

x  y. 

o 


(245c) 

(245d) 
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Ayyyy  _  2(A  y  +  6A  y  y} . 

o  o  o 


(245e) 


It  is  worth  noticing  that  Eqs.  (241) -(244)  are  valid  for  schemes 

t  t 


of  higher  degree 
t  t 

in  the  operators  t  ,t  ,d  ,  and  d  ,  as  long  as 
x  y  x  y 

X  X 

A 

=  A  y  y  =  0,  i.e. 

,  as  long  as  composite  transport  is  excluded. 

With  the 

above  equations, 

Eqs.  (231)- (233)  yield 

t  t 

(log  A)  =  iA  ,  (log  A)y  =  iA  y ; 

o  o  o  o 

(246) 

(log  A)XX  =  - 

d  t  _ 

■  2A  X  +  (A  X)2; 
o  o 

(247a) 

(log  A)Xy  =  - 

t  t  t  t 

■  A  X  y  +  A  XA  Y; 
o  o  o 

(247b) 

(log  A)^  =  - 

a  t 

■  2A  Y  +  (A  Y) 
o  o 

(247c) 

(log  A)XXX  = 

t  d  t  d  t 

x.  .  ...  x  x  .  x.  3. 

-  lA  (1  -  6A  )  -  x(6A  +  2 (A  )  ] ; 

o  o  o  o 

(248a) 

(log  A)XXy  = 

dt  dt  ttt  tt 

x  y  „  x,  y  x,.xv  ,  x,  y. 

-2i(A  ^  -  A  A  ^  +  2iA  (A  *  -  A  A  ) ; 

O  OO  o  Q  oo 

(248b) 

xyy 

(log  A)q 

dt  dt  ttt  tt 

-  2i  (A  y  x  -  A  yA  x)  +  2iA  y  (A  X  y  -  A  X.A  y)  ; 
o  o  o  o  o  o  o 

(24P  '.) 

(log  A)yYY  = 
o 

t  d  t  d  t 

-  iA  y (1  -  6A  y)  -  i [6A  y  y  +  2 (A  y) 3] . 
o  o  o  o 

(248d) 

Only  the  odd  derivatives  are  imaginary.  Therefore,  Eq.  (229)  implies 


(log  A) ' 


(log  A)y  (log  A)XXX  S3 

- — 2.  B  ]  +  [ - ; — 2_  -Ji 

i  y  i  6 


do,  A|xyy 
o  x  y  o 


5  S2  (log  A)yyy  B3 

*  y  +  _ _ _z 

2  i  6 


- —  ]  +  ...  (249) 

6 


Example : 


Let  us  analyze  the  phase  error  associated  with  Eq.  (226e)  , 


A  =  (1  -  e  t  -  e  t  )  +  v  d  +  vd  -  pdd-et)-ud(l-et) 
xx  yy  xx  y  y  x  x  xx  yy  yy 
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1 


By  direct  differentiation  we  get 


t  t 

,  x  .  y 

A  =  £  ,  A  =  -  e  ; 
oxo  y 

d  d 

.  x  .  y 

A  =  v  -  y  ,  A  =  v  -  y  ; 
o  x  x  o  y  y 


A  =  0,  A 


£  y  i 

XX 


t  t  y  y 

.  y  y  „  A  =  £  y  , 
A/  *  =  0,  o  ypy 


d  d 

A  X  X  =  0; 
o 

d  d 

A  ”  =  0; 


tt  td  td  dd 

A  X  *  =  0,  A  X  y  =  0,  A  y  X  =  0,  A  X  * 
o  o  o  o 


whence 


(log  A)' 


(log  A) 


(log  A)' 


(log  A)- 


=  -  V 


=  6ex<!-  vx  +  f) 


=  -  2e  (v  -  y  -  £  ) ; 
y  x  x  x 


(log  A)XX* 

O  _  .  2 

- : -  =  -  2c  (v  -  y  -  C  )  ; 

l  x  y  y  y 


(log  A)fy  z2 

- : — 2 —  =  6e  (-  -  v  +  -?)  . 

i  y  6  y  3 


Substituting  in  Eq.  (249) ,  we  get 


6  "  exBx  [1  +  (v  -  -  -  ^)82  +  (y  -  y  -  £2)62  + 

x  6  3  x  y  y  y  y 


e  8  [1  +  (v  -  ~  -  -g-  )8  2  +  (v  -  y  -  e  2)S  + 

y  y  y  6  3  y  x  x  x  x 


(252a) 


(252b) 


(252c) 


(252d) 


•vc  * 

0,  •** 
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Noting  that  e  8  =  (—7 — )  (k,  6x)  =  u,k.  c't  and  e  6  =  u„k_6t,  we  can  rewrite 

xx  oxl  11  y  y  2  2 


the  above  equation  as 


d  =  -  k  •  V<$t, 


where  k  *  (k  ,k0) ,  V  =  (v^,v2>.  If  U  =  (.u^,u^) , 

_2 

V  =  u,  [1  +  (v  -  \  -  —)  &2  +  (v  -  U  -  £2)  8^  +  . 
i.  1  X  o  3  X  y  X  y  y 


v2  =  u2  [1  +  (vy  -  I  -  (VX  "  px  ■  £x)Bx  +  •• 


(254b) 


Ccomparing  Eq.  (253)  to  the  analytical  solution,  we  find 


d  .  .  =  -  k  •  U6t 

analytic 


as  is  obvious  from  Eqs.  (211) . 


Following  Eq.  (19b),  we  define  a  relative  phase  error  matrix,  R, 


such  that 


V  =  U  +  RU, 


where  R,  given  by  Eqs.  (254),  in  this  scheme  is 


(v  -  \  -  e2/3)62  +  (v  -  y  -  £2)32  ... 
x  6  x  x  y  y  y  y 


1  2  2  2  2 
(v  -  --  e/3)B-Kv  -  y  -  £  )  8  +  . 

y  6  y  yx  x  xx 


Thus  we  can  reduce  the  phase  error  to  fourth  order  by  selecting 


1  2  ,, 
v  =  —  +  £  /  3 ; 
x  6  x 

y  y 


(258a) 


v  -  u  =  £ 

XXX 

y  Y  Y 


(258b) 
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Solving  Eqs.  (258a,  b) ,  we  get 


12  2 
U  =  T  ~  ~  £ 
x  6  3  x 

y  y 


(258c) 


AMPLITUDE  ANALYSIS 


Following  the  analysis  of  the  one-dimensional  case,  since  t  A I  =  AA 


we  have 


(|aD0  =  1» 


( |a| 2)X  =  AXA~  +  A  AX 
o  o  o  o  o 


( I A I 2) y  =  AyA  +  A  AY, 
o  o  o  o  o 


and  so  on.  Noticing  from  Eqs.  (241M245)  that  odd  derivatives  are  purely 
imaginary  while  even  ones  are  real,  we  get  after  substituting  in  Eqs.  (259) 

( | A | 2 ) o  =  1;  (260 


(|A|2)XX  =  2(AXX  +  (f  )2]; 


AX  AY 

2[Axy  +  -4  -4i, 

O  11 


^|2>oY 


(261a) 


(261b) 


(|A|Vy  -  a  [A™  .  t-r-)2] , 

O  O  1 

Axxy 

(|A|VX**  .  2[A™“  .  4  (“  )  (—7 - )  +  3  (A**)  2  ]  ; 

O  O  11  o 


„ xxx  . y 
A  h 


AXXY  AX 


(|A|2)XXXy  =  2(AXXXy  +  (^~)(4>  +  3  (— — — )  (4-)  +  (3AXXAXy); 


(261c) 


(262a) 


(262b) 


AXXY  Ay 


(  !  A  j  2 )  ^XyY  =  2  [AXXyY  +  4~)4)  +  AXXAyy 


AXyy  AX 

+  2  (— -r  )  (-r-  )  +  2(AXy)2]; 


(262c) 
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<iAr> 


2,xyyy 


A*  Ayyy  A*yy  Ay 

2  [AXyyy  +  (-£)(-£r-)  +  3  )  +  SA^A^l, 

O  XI  XI  O  O 


(lAT) 


2,yyyy 


Ay  ayyy 

2[Ayyyy  +  4  (— r-  )  (— r - )  +  3(Ayy)2], 

o  ii  o 


while  the  odd  derivative  vanishes.  Consequently,  the  expansion  of  j  ? 
takes  the  form 


wv  +  (i'i2>:yVy  + 


S3e 


2  2 
8  3 


<•  ,w2)ry  dAi2)^^ 


*  ,|a]2>™^+  * ... 


Example : 


Using  Eqs.  (250),  Eqs .  (241H243)  yield 


Aq  -  1; 


X  y 

A  A* 

Q  _  __0 

i  x '  i 


*xx  «•>/  x 

A  -  -2(v  -  U  ) 

O  XX 


1 

-2v 


AXY  =  0; 
o 


AYY=  -2 (v  -  u  )  = 
o  y  y 


'2V 


which  when  substituted  into  Eqs.  (260) ,  (261)  give 


<l*r>0  = 


[  A I  2 ) xx 

1  1  o 


2e2; 

x 


3Ba 


(262d) 

(262e) 


(263) 


■  2c  £  ; 


dsl2)"  .  -  2c2. 

1  1  o  y 


Equation  (263)  then  expresses  the  amplitude  expansion  as 


|a|2  =  1  -  (-V  -  2e  e  g  3  +  e2B2)  +  ... 

1  1  xx  x  y  x  y  y  y 


|a|2  =  1  =  -  |  e  8  -  e  6  |2  +  . . . , 

ii  1  x  x  y  y 


showing  that  the  diffusion  of  the  scheme  is  second-order. 


POSITIVITY 


From  Eqs .  (223)  and  using  Eq.  (258a) 

»I?j  -““.i11-  2(i  +  T'  -  2(l*f>' 


n  r  .1  ,  '"‘x  ,  lx,  ^  n  rl  .  -x,  ^  Ex  . 

+  °i+l,jt(6  +  T>  -  T1  *  6  “3  T1 


-  f'  'f' 


e  e 


Each  of  the  square  brackets  is  >_  0  for  |ex|,|e  1^ 

TD  n  „ 

Consequently,  ^  >_0  if  all  ^  0 .  Now  we  get 

ijt  T  T  T  T  T 

n+i  td  .  x  x  t  x  ,  .  ,  y  -  y  ^  „  y  > . 

p.  .  =  p.  .  -  u  (p.,-1  .  -  2p.  .  +  c.  ,  .)*  -  y  (p.  .  ,  -  2p.  .  +  p.  .  ,) 

1,3  1,3  x  1+1,3  i,3  1-1,3  Y  1,3+1  1*3  i,3“l 

The  asterisks  denote  the  fact  that  the  antidiffusion  fluxes  are  trimmed 

enough  such  that  pn+^  is  limited  by  the  sign  of  Then  pn+^  >  0. 

i,j  1,3  i  ,  j  — 


f  "  * — f  x  c 
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STABILITY 

Equation  (264)  proves  the  stability  of  the  scheme  near  6^ 

For  the  scheme  to  be  completely  stable,  however,  we  must  have  ]a| 

0  <  3  , 6  <  7T.  Let  us  check  A  ,  at  the  largest  values  admitted 

-  x  y  - 


namely 

1/2. 

From  Eq.  (226e) 

ar 

=  1 

-  2(vx  -  „X>(1 

-  cos  6  )  -  2  (v 
x  y 

u  ) (1  -  cos  6  )  ; 

y  y 

AI 

=  - 

e  s in  6  [  1  + 

X  X 

2u  (1  -  cos  6  ) ] 

X  X 

- 

£  sin  6  [1  + 

y  y 

2u  (1  -  cos  6  )  ] 

y  y 

Subs ti tut in 

.g  for  v  -  y 

XX 

y  y 

from  Eq.  (258b) ,  and 

U  from  Eq.  (258c) 
X 

y 

ar 

=  1 

0 

-  2e“  (1  -  cos 

X 

Bx)  ~  2e^(l  -  cos  6v 

)  ; 

AI 

=  - 

£  sin  6  [1  + 

X  X 

4(1  -  4e") (1  -  cos 

4  x 

6x)] 

- 

£  sin  6  [1  + 

y  y 

4(1  -  4e“) (1  -  cos 

4  y 

V1* 

At  £  , 

X 

“y 

'  1/2,  Ux  =  Uy 

=  0;  Eqs.  (267)  reduce  to 

ar 

=  i 

1  -1  J 

-  —  a  -  cos  6 

2  x 

)  -  4d  “  c°s  ; 

ai 

=  - 

~(sin  6  +  sin 

X 

v- 

2 

Noticing  that  1  -  cos  8  =  2  sin  (6/2)  and  sin  6=2  sin  (6/2)  cos 


r- 


*6*0. 

y 

<  1  for 

for  £  ,  £  , 

X  y 

(266a) 

(266b) 

we  get 

(267a) 

(267b) 


(6/2)  , 


we  get 


72 


A  =  1  -  sin2 (3  /2)  sin2  (3  /2); 
R  x  y 


A  =  -  [sin  (3  /2)  cos  (3  /2)  +  sin  (3/2)  cos  (3/2)), 

x  xx  y  y 


yielding 


.  . 2  2  2  .  4  6x  .  4  ^  .  2  Bx  .  2  j 

|  A  |  =  A  +  A  =  1  +  sin  —  +  sin  +  2sin  —  sin 

K  I  2.  z  4  2 


2  Bx  .  2  Sy  .  2  Bx  2  Bx  .  2 

-  2  (sin  —  +  sin  )  +  sin  —  cos  —  +  sm 


cos 


2  !z 


x  y  x  y 

+  2  sin  —  sin  cos  —  cos  . 

2  Bx 

Collecting  the  terms  containing  sin  —  ,  we  have 

-  sin2  -y-[(l  -  sin2  -1 )  +  (1  -  cos2  —■ )  }  =  -  sin2  ytcos2 

.  2  Bx  .  2  Bx  2  XL  •  4  8x 

+  sm  —  J  =  -  sin  —  cos  -  sm  — . 


6  3  3  3 

2  y  2  y  2  x  4  y 

Similarly  sin  terms  yield  -  sin  -j-  cos  —  -  sin  —  ,  resulting 


in 


i  |2  2  Bx  2  By  .  Bx  .  Bx 

| A |  =  1  -  (sm  —  cos  2  ~  2  Sin  ~  Sln  2  cos  ~  cos  , 


S 

JL 


,3  _  S 

•  2  y  2  x  . 

+  sm  n*-  cos  — )  = 


6  3 

x  y 

-  mQ  — ±-  — 

2 


S 


x  ,  2 


i  •  2  3x  2  Bx 

-  1  -  sm  (—  — ■£■)  =  cos  (—  — ) . 


Consequently, 


ex=ey=i 


,  x 

cos  (— 


e .  e 


(268) 


showing  the  scheme  to  be  stable  at  e^,  =  J.  The  value  J  A |  for  smaller 

values  of  was  evalu  ted  numerically  and  found  always  to  satisfy  <_  1. 

Hence  the  scheme  is  completely  stable. 


4 
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It 

In  fact,  on 

to  the  3  = 

x 


is  worth  noticing  the  diagonal  symmetry  of  Eas.  (264) ,  (268) 


the  3^,  plane,  |Aj  looks  like  a  wave  with  front  parallel 
6^  diagonal,  as  illustrated  in  Fig.  27. 
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XIII.  RECTANGULAR  GRID  MOTION 


Consider  a  system  of  points  tagged  by  the  double  indices  i,  j ; 


x.  =  x (i, j ,t) ; 
1  >  j 


(270a) 


yi(j  -  y(i»j/t) 


(270b) 


Figure  28  illustrates  the  grid  formed  by  Eg.  (270)  at  a  given  time  t. 
The  pair  of  numbers  at  each  point  indicates  (i, j) .  For  a  strictly  rectangular 
grid  at  all  times  (which  includes  Langrangian  grid  motion) , 


x  -  x(i,t) ; 

If  J 


(271a) 


y,  ,  =  y(j/t) . 

1  >  J 


(271b) 


which  we  therefore  denote  from  now  on  by 


x  =  x  (i,t) ; 


(272a) 


y  j  =  y(jft) . 


(272b) 


This  leads  to  a  mesh  as  in  Fig.  29. 


*•»  ^  'n 


v-  - 
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XIV.  GEOMETRICAL  ASPECTS 


We  consider  seven  geometries.  (These  by  no  means  cover  the  whole 
spectrum  of  two-dimensional  systems.)  In  cartesian  geometry,  we  have  x-y 


(x-z  or  y-z) ;  in  cylindrical  geometry,  r-z,  r- p ,  and  z-j,  and  finally  in 
spherical  geometry  r-c ,  r-p,  and  5-?.  Figure  30  illustrates  a  finite 
control  volume  in  each  of  the  different  cases. 

As  explained  earlier,  when  the  grid  moves  the  control  surface  area 
in  the  integral  form  of  the  conservation  equations  should  be  an  average 
surface  area  defined  as 


(u  '5t#n)dS  =  swept  volume. 


In  one-dimensional  cases,  the  above  definition  reduces  to  defining 


A . 

interface 


swept  volume 


.  In  two  dimensions,  however,  this  is  not  enough. 


We  have  to  find  a  path  between  the  old  grid  and  the  new  one  such  that  we  can 
construct  a  mean  ceil  having  its  surfaces  equal  to  the  average  areas  and 
corners  located  on  the  above  path. 


1.  Cartesian  Coordinates 

Figure  31  illustrates  the  location  of  cell  (i,j)  at  times  tn  and 
tn+1  =  tn  +  6t.  The  left  and  right  interfaces  are  denoted  by  (i  -  1/2,  j) , 
(i,j  +  1/2,  j),  respectively,  and  the  bottom  and  top  ones  by  (i,j  -  1/2), 


(i,j  +  1/2).  We  notice  here  that  since  all  i  ±  1/2,  j  interfaces 

(different  j's)  move  as  a  whole,  the  grid  velocity  is  independent  of  j. 

It  is  therefore  denoted  by  u  .  without  a  j  index.  The  same  is  true 

i±i 

for  v?  . , . 

3±i 
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In  cartesian  geometry,  it  is  obvious  that  the  path  needed  is  a 
straight  line  between  the  new  and  old  corners  of  the  cell,  and  the  mean  cell 
is  halfway  between  the  old  cell  and  the  new  one. 

The  volume  swept  by  interface  (i  ±  1/2,  j)  is  given  by  the  product 
(average  base)  x  (height) : 

,  n  n  .  ,  n+1  n+i. 

lij.j  2  111  iil  *iii' 

1 /„n  ^  .n+1. 2  n+1  n  , 

=  —(A.  +  A.  )  x.  ,  -  x.  , 

2  j  j  i±i  i±i 

where  we  notice  again  that  the  i  index  is  omitted  from  A.  ,  .  since  all  the 

i±i ,  3 

A^+j  interfaces  (different  i's)  are  equal.  The  above  equation  can  be 
written  as 

n  n+1  n  n+1 

showing  that  the  mean  area  is  halfway  between  old  and  new.  The  grid  velocity 

g 

u^+j  in  this  case  is  considered  constant, 
n+1  n 

„g  _  *i±L " 


and  the  mean  area  is 

An+i  =  n+i  _  n+i 
D  yj+i  yj~i 


.n  .n+1 
A  .  +  A  . 


(288a) 


where 


n+i  1 ,  n  n+1. 

yj±i  =  2<*j±i  +  ^j±i) 


(288b) 


Similarly, 


„n  n+1 

,  .  ,  A .  +  A . 

n+i  n+i  n+i  i  i 

A.  =  x.  .  -  x.  ,  =  - - - 

1  i+i  i-i  2 


(289a) 


where 


n+i  1,  n  ,  n+1 

Xi±i  2(xi±i  i±i}‘ 


(289b) 


The  mean  cell  colume  is 


n+i  n+i  n+i  n+i  n+i 

V.  .  =  (y  .  ,  -  y  .  , )  (x .  .  -  x .  . ) 

3+i  3~t  i+i  i-i 


(290) 


2.  Cylindrical  (r-z)  Coordinates 

Let  us  derive  the  required  path  between  the  corners  of  old  and  new 
cells  such  that  the  corners  of  the  mean  cell  fall  on  that  path.  Figure  32 
illustrates  the  old  and  new  cells.  Figure  33  shows  the  volume  swept  by 
interfact  (i,j  ±  1/2): 


n+1 

2j±i 


¥. 


i,j±i 


*rLi  dh±i  -/ 


n+1 

Zj±i 


2 

irr.  ,  dz..,, 
i-i  3±i 


'jii 


'j±i 


(291) 


where  it  is  obvious  that  a  linear  average  can  be  obtained  if  r.,,  is 

i±i 


assumed  to  be  linear  in  z .  . .  Let 

3±* 


,  n+1, 2  ,  n  ,2 

"it!*  -  ‘Lit1 


,  n  ,  2  n 

ri±i  =  zjti  "  z j±i 

n+1  n 
Z3±i  •  Z3±* 


(292a) 


i.e.,  a  parabolic  path.  The  above  formula  can  be  written  concisely  as 
2 


A  r 


i±i 


dz  . 


“Li 


i±i 


(2 
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yielding 

r2  =  (rn  )2  ^  AR2 

i±i  1  i±r  '  AZj±i  i±i 
Substituting  (292c)  into  (291) ,  we  get 

r1  2  2  ^Zi  +  i 

AVi,j±i  =  "AZj±i  1  (ri+i  ‘  ri-i)d(Ii^° 


.  n  ,  ,  n+1.  2  ,  n  .2  n+1  2 

,lri+J  +  !r^i>  ,  l*W  ♦  “W  , 

■  *“j*l 1 - 2 -  +  - 5 - 1  ' 


yielding  a  mean  area 


A» 

An+*  =  =  tt  [  (rn*t)  2 

1  AZj±i  1+i 


.  ,  An  +  An+1 
(ri-i5  1  -  2 


where 


,  n  .  2  ,  .  n  2\i 

(ri±ii  +  (ri±i} 


(292c) 


(293a) 


(293b) 


This  shows  the  advantage  of  the  parabolic  path  (292a)  ,  namely 


Azni  atiH 

AZ  2 

jti  AE.ti 


=  1/2, 


(293c) 


i.e.,  the  average  area  A^  is  halfway  along  z  between  the  old  and  new  ones, 


n  ,  n+1 

n+i  Zj±i  2j±i 


(293d) 


Following  the  nomenclature  of  Fig.  34,  the  volume  swept  by  the  inter¬ 


face  i±J,j  is 
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A*. 


n+1 

r 

'i±i 


iti.j 


2"ri±i(zj+i  “  zj-i)dri±i 


'i±i 


n+1 
'  riti 


n 

ri±i 


(294) 


But 


j±i  Zj±i  +  Azj±i' 


(295) 


yielding 


Av. 


i±i » j 


n+1 

ri±i 


'i+i 


Az .  .  Az  .  , 

wl(lj+4  ~  +  AVi  “  AVi)ld<±i 


which,  with  (292b)  gives 


AV. 


i±i,j 


jD2  r ,  n  n  , 

nA\±i  [U3  +  i  ~  Vi* 


0 


Ar2  Ar2 

+  '^-i  (AV* ' 


AR 


i+i 


AR' 


i±i 


(296) 


Now  to  be  able  to  construct  a  rectangular  mean  cell  with 
the  parabolic  paths  of  Eq.  (292a) ,  the  quantities  j 
half  way  between  old  and  new,  i.e.. 


its  four  corners  on 
have  to  be  also 
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where  we  notice  that  the  interface  area  is  dependent  on  both  i , j ,  in  contrast 
to  the  cartesian  case.  With  Eq.  (296) ,  Eq.  (297)  yields 


A¥i±i 


i±i,j  AR2  /2rn+i 
i±i/2ri±i 


whence,  if  denotes  the  average  velocity  of  the  grid  during  t , 

n+1.2  .  n  .2 

«  ■  'l&Lip* iL 


where,  as  is  clear  from  Eq.  (293a) ,  it  was  assumed  that 


9  ,  n+1  n 

vjtiSt  *  4zj±l  '  zj±i  -  2jii' 


The  difference  between  the  form  of  (299)  and  that  of  (300)  is  attributed  to 
the  parabolic  path  of  the  corner.  If  the  grid  velocity  vg+^  is  a  constant 
during  fit,  we  evaluate  ug+^  at  tn+^  =  tn  +  from  (292c) 

dr  ,  ART  ,  dAz  ,  AR2  ,  dz  , 

-  rti  i±i  _ l±i  x±i  j±j 


"i±i  dt  AZ  dt  AZ 

3±i  3±i 


dt  ' 


whence 


AR2 

i±l 

vg 

j±i  . 

ARi±i 

2ri±i 

A2j±i 

'  2ri±i 

using  Eq.  (300).  Since  v...  =  const.,  (A z. . , /AZ . ) 

3±i  J1:  t  ot 

Consequently,  from  (292c) 


2  |  lr,  n  .2  ,  ,  n+1.2.  ,  n+i,2 

lit1,,  it’  2I(W  *  “W  1  ■  lri±l>  ’ 
2 


=  1/2. 


81 


thus  reproducing  Eq.  (299)  when  substituted  in  Eq.  (301).  A  more  general 
definition  of  the  average  interface  is  therefore 

[" mean  1  swept  volume 


interface  (velocity  of  interface  at  tn  +  ^p).5t 
area 


where  the  denominator  is  approximately  but  not  quite  exactly  equal  to  the 
distance  the  interface  is  shifted.  Finally,  the  mean  cell  volume  is 


-  „f,  n+i,2  ,  n+}  2.  .  n+i  „n+i 

vi,j  -  ,I0W  '  °W  <‘v»  - 


3.  Spherical  r-c  Coordinates 

Figure  35  illustrates  cell  i,j  at  tn  and  tn+1.  Consider  Fig.  36 
showing  the  motion  of  interface  i,j±i.  The  volume  swept  by  interface 
i.j±i  is 


AV.  ... 
i/3±i 


2ti-  ,  3  3  . 

T  (ri+i  -  W  sin  0j±i  d0j±i 


2  n  3  3  ,  .  .  .  , 

T  (ri+i  "  ri-i)Q(-  COS 


showing  that  we  can  get  a  linear  average  if  r^  +  ^  is  assumed  to  be  a  linear 


function  of  (-  cos  9^  +  ^)  •  Let 


(w 

3 

/  n 

3 

-  n 

(JW 

cos  6  . 
3 

n+1 

'he1 

3 

"3 

(ri±j} 

cos  e . 

r 

,n+l  ' 


(305a) 


*  *  “*  *■  “  A*-  *• 
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or  in  a  more  concise  form. 


Ari±i  _  A(-  cos  e.±i) 
AR3+i  ~  A<" cos  ‘W 


(305b) 


This  yields 


ri±i 


„  ,  A  (-  cos  0  .  ,, ) 

(rn  )  3  +  _ 1*L_ 

i±i'  A(-  cos  ®j+^ 


aEiir 


(305c) 


Substituting  into  Eq.  (304)  we  get 

-f4'-  =°s  wf  UrlU  '  ri-i’ 


A(-  cos  0  ,) 

3  ~  *  /a  r>-^ 


+  — - (AR3  -  AR3  ,)]d 


A(-  cos  6^  +  ^) 


A(-  cos  0.  ,)  vtli+4  i-4'J“  A  (-  cos  0.., 

3*5  3*4 


2  TT 


n+1. 


=  —  (cos  5  ^  -  cos  9  ^)  t 


,  n  3  .  ,  n+1  3 

(1W  * 


,  n  ,  3  ,  .  n+1.  3 

<ri-j’  *  <ri-i’  , 


(306) 


We  notice  that  the  mean  interface  i,j±4  is  halfway  on  a  cosine  scale  between 

„n  n+1  ,  .  ,  n  n+1  „ 

0  »  9  .^.1  or  on  a  cubic  scale  between  r .  ,  r.^,.  As  for  interface 

3±4  3±4  i±4  i±4 

(i±4,j),  it  sweeps  a  volume  (see  Fig.  37), 


n+1 

ri±i 


AV. 


i±4/ j 


n 

ri±i 


'  n 

ri±4 


27rr.,,(cos  0.  ,  -  cos  e.^.Jdr... 

i±4  3-4  3+4  i±4 


n+1 

ri±4 


j(cos  -  cos  0j+i)dr3±r 


(307) 
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But  cos  9.  .  =  cos  6 ,  .  +  A (-  cos  9  .  ,  . ) ,  yielding 
3±i  3±i  3*4 

n+1 


A(-  cos  9.+  ,  ) 

[(cos  9n  -  cos  an  )  +  — — — { A(-  cos  0.  .) 

3“i  3+4  £(-  cos  0j+^)  3-1 


i±i,j  3 


-  A(-cos  0.+^)}]  dri  +  ^' 


which  with  Eq.  (305b)  results  in 


AV.  ,  .  =  —  AR3  , 

i±4,j  3  i±4 


,  ,  .n  „n+l, 

[ (cos  6 .  ,  -  cos  3  .  , ) 

3-1  3+4 


,  3  3 

Ari+i  Ari  +  1 

+  - {A<-  cos  e.  ,)  -  A  (-  cos  0.  ,)}]  d  — r-^- 

3'*  3  * 

n  .n+1  n  .n+1 

_  .  cos  0  .  ,  +  cos  9 .  ,  cos  6 .  ,  +  cos  9 .  . 

«  ar3  [  _ H* _ Hi  -  _ a±i _ Hli 

3  i±4  1  2  2 


Here  we  notice  that  interface  i±i,j  is  halfway  on  a  cubic  scale  between 

n  n+1  ,  n  .n+1 

r .  r  .  J  .  or  on  a  cosine  scale  between  9..,,  9.  Consequentlv ,  we  can 

i±4  i±l  3±4  3±4 

construct  a  mean  cell  having  its  corners  on  the  paths  of  Eq.  (305a) .  Let 

/  n  ,3  n+1  3\  1/3 

.... ,  Km!  ; 


,n  „n+l 

,  i  cos  0  .  ,  +  cos  9  .  , 

.n+4  _  j±i  1+4 


Eqs.  (306)  and  (308)  can  be  written  then  as 


.  2  it  n  .n+1.  r  ,  n+4>  3  ,  n+i 

i,  j±4  =  T(cos  ej±4  ~  COS  3j±4)  ri+4  ‘  (riH 


and 


(312) 
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2tt  n+i  ,n+i  n+1  3  ,  n  ,3. 

AVi±i,j  =  T(cos  ej-i  -  cos  ej+i)C(ri±i)  "  (ri ±i>  ]' 

respectively.  Now  in  order  to  be  able  to  construct  the  average  cell, 

A.  .  .  and  A .  .  .  .  should  take  the  forms 
1,3+i  i±i,3 

„n+i  .  .  n+i. 2  .  n+i. 2,  ,n+i 

Ai,j±j  '  *l(1W  -  1  sln  ”j.r 


forcing  the  choice 

Av 


vg  fit  =  . 

i,l±i  ,n+i 

Ai,j±i 


.  n+i. 3  ,  n+i  3  .n  _n+l 

'u  o  (r  ,)  ”  (r.  .)  cos  a...  -  cos  9  .  ,  , 

iQ±i  _  2  l+i _ iyi  3±i _ T±i 


3  ,  n+i. 2  .  n+i,  2  n+i 

lr1+t>  -  "i-j1  Sln 


and 


n+i  _  .  n+i, 2,  „n+i  „n+i, 

A.,,  .  =  2i(r.  ;)  (cos  9.  ,  -  cos  0  ?) , 
i±i»3  i±i  J-i  3+i 


forcing  the  choice 

A*i*i,i  =  (giti>- 


n+l,  3 


Ui±*  5t  = 


An+i 

i±i»  j 


3(r"t»)2 


To  complete  the  formulation,  it  remains  to  check  the  consistency  of  the 
g  g 

velocities  u^+^,  vj_  +  namely  that  they  occur  at  the  same  instant. 

Differentiating  (305a)  with  respect  to  time  and  taking  r.+,,  6,+,  at  the 

.  ,  .  n+i  Qn+i  1_  3~ 

moment  wnen  they  are  halfway  i.e.,  r^  +  ^,  9j  we  get 

3  (ri±j )  "  (dri±^/dt)  n+*  _  sin  9*+j  (d9 j±i/dt)n+* 

(ri±i1)3  -  (<+i)3  c°s  e*ti  -  cos 


Recognizing  that  (dr^  + ^/dt) =  ug+^,  the  velocity  of  the  grid  at 

.  n  <St 


t  + 


2  ' 


from  Eq.  (314b)  and  the  above  equation  we  obtain 


d9  . 


n+i 


j±i 


dt 


-)•  St  = 


cos  9 n  .  -  cos  9n"*"i 

_ l±i 

Qn+i 

sm  9  ,  , 

:±i 


n  6t 

whence  from  Eq.  (313b) ,  the  velocity  of  the  grid  at  t  +  —  , 


(313a) 


(313b) 


(314a) 


(314b) 

two 


(315) 
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„  ,  (r^y  -  <rny  d5  . 

1/jii  3  .n+i  2  n+i  2  dt 

1  i+r  “  1  i-i' 

dej±i 

Recognizing  ( — 4- — )as  tne  angular  velocity  at  t  t  ot/2  of  the  interface 
at 

i,j±i  [from  Eq.  (315)  obviously  independent  of  index  i  as  expected  after 
the  discussion  in  the  section  "Rectangular  Grid  Motion’’],  we  can  define 
an  average  radius  for  the  interface  i,jii  (independent  of  j)  as 


*  <!>> 


Finally,  the  mean  cell  volume  is 


n+i  2 tr  n+i  3  ,  n+i,  3.,  „n+i  .n+i, 

h,3  =  TI(1W  -  1(cos  V*  -  “s  Vi* 


Coordinate  cases  4-7  will  be  treated  in  a  later  report. 


W*  * 
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XV.  SOURCE  TERMS 


As  explained  earlier,  source  terms  are  integrated  either  over  cell 
volume  or  over  cell  interface  area.  The  volumes  and  areas  used  are  those 
of  the  mean  cell.  The  balance  of  source  terms  is  the  main  reason  for  the 
necessity  of  a  closed  mean  cell  construction,  i.e.,  the  ability  to  construct 
a  closed  cell  whose  corners  are  on  the  paths  between  old  and  new  cell 
corners.  For  example,  if  we  try  to  solve  a  hydrostatic  pressure  problem 
in  cylindrical  r-z  coordinates,  the  momentum  equation  is  nothing  but  the 
balance  of  the  body  gravity  force  and  the  pressure  force  on  the  top  and 
bottom  surfaces.  If  the  mid-cell  interfaces  compose  a  closed  surface 
enclosing  the  midway  cell  which  happens  to  have  a  volume  consistent  with 
the  interface  areas,  force  balance  is  already  guaranteed  (provided  the 
pressures  are  correct) . 

Let  us  consider  the  difference  form  of  -  I  pnds  (yielding  -  grad  p) 
in  tne  three  coordina*  ises  considered  above.  The  resulting  forces 
along  the  x,  y  directions  are 


F 


x .  . 


,  n+i 


P 


n+i 
i+i  /  j 


(P 


n+i 

i* 


P 


n+i 
if  j+i 


respectively,  where 


n+i  4.  r,n+* 

n+i  -  Pi.j  Pi+lfj 

Pi+i,j  2 


(321a) 


(321b) 


(322a) 


and 


i»  j+i 


n+i  n+i 
Pi,i  +Pi,1+1 


(322b) 


In  cylindrical  r-z  coordinates, 

_  n+i  ,n+i  n+i  n+i  _  n+i.  n+i  n+i,  ,  n+i  n+i 

F  =  P .  :  .A.  I  .  -  p.  i  .A.  1  .  +  2-  p.  r.  ,  -  r.  Hz.  -  z .  , 

r  .  i-i,3  t~i.3  i+i.3  i+i,3  a ,  3  i+i  i-i  3+i  ]-« 

1  •  3 


which  with  Eq.  (303)  yields 


pn+i  ¥n+i 


F  =  pn+J  .  An+J  .  -  p n+\  .  An+\  .+  1^ 

r.  .  l-i , 3  i-i/3  i+i. 3  i+i #3  n+i 


(323a) 


n+i  n+i 
i  r.  -  +  r.  : 
n+i  _  i+i  r-i 


We  notice  that  p^+^/r^+^  acts  as  a  body  force  per  unit  volume.  The  force 
in  the  z-direction  is 


_  .  n+i  n+i  ..n+i 

F  =  (p.  .  -  -  p.  .  , )  A . 

2.  •  i » 3“i  i.3+i  i 

i  '  3 


(323b) 


In  spherical  r-6  coordinates,  as  illustrated  in  Fig.  38,  the  pressure 
acting  on  the  hatched  area  creates  a  resultant  force  normal  to  the  axis 
from  which  3  is  measured,  which  in  turn  gives  rise  to  a  radial  component 

I  I 

F^  and  a  tangential  component  F„.  This  situation,  namely,  the  creation 
of  a  body-force- like  component,  occurs  whenever  the  area  of  parallel 
surfaces  of  the  cell  are  not  equal.  This  is  bound  to  happen  whenever  the 
interface  area  depends  on  both  indices  i,j.  A  simple  way  to  evaluate 
the  force  generated  is  the  "pressure  x  projected  area"  since  this  area 
is  the  difference  between  the  areas  of  these  parallel  surfaces.  The 
radial  force  is  therefore 


cm***- 


.  A 
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F  n+i  n+i  n+i  An+i  +  F* 

ri,j  '  Pi-i,j  Vi,j  Pi+i/3  i+».3  ri,j 


where 


it  j 


n+i. 2, 


„n+i 


2pn+i  V?+i 
r  i  ^  ■»  .  1 


n+i.  _  1 iq  J-rJ-  . 


(r^D  (cos  e".j  -  =os  0j+i)  -  —-n+i 


R.  was 
1 


defined  earlier  in  Eg.  (317).  Thus 


2  n+i  ¥n+i 

n+i  an+i  _  n+i  An+i  _  +  l.j  . 


Fr  =Pi-i,j  i'i-3  Pi+i.3  i+i<3 
i/3 


.n+i 


As  for  the  tangential  direction, 

n+i  An+i  -  Pn+i  ,  An+^,  +  P* 

■ i,  j-i  i/j+i  t/j+i  1»3  i  i,j 


F .  =  P, 

i/ j 


where 


n+i  ,.n+i  _  An+i  ,  .  ,  pn+i  [(rn+|>2 

=  Pi,j  Ai. j+i  i/3'i  Pi'3 


3 .  . 

i/l 


-  (rj**)2]  (sin  ~  sin  5 • 


If  we  note  that 


n+i  an+i 

sin  ej+i  -  sm  e._i 


=  2  sin 


an+i  _  .n+i 
( h±Ll—ti)  cos  (■ 


n+i  -n+i 

iti  Hid) 

2 


and 


,n+i  +  0n+i 


cos  0^  -  COS  -  2  sin  (- 


,n+i  .n+i 
2 


and  using  Eq.  (318) ,  F  can  be  expressed  as 
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XVI.  ALGORITHM 


We  describe  the  implementation  of  the  scheme  of  Egs.  (223) .  The 
program  and  calling  sequence  are  listed  in  the  Appendix. 

Assume  a  rectangular  grid  in  two  dimensions  denoted  by  the  coordi¬ 
nates  x,  y  (not  necessarily  cartesian;  for  example  x  E  r,  y  =  6  yield 
spherical  coordinates).  Let  the  interfaces  coordinates  be  X3/2'‘‘*' 

\+l/2;  yi/2'  y3/2 . yNy+l/2  (S6e  Fig*  39)- 

The  cell  centers  are  located  midway  between  the  interfaces  and  are 
denoted  by  a  pair  of  indices  (i,j),  corresponding  to  (x,y),  respectively. 
The  cell  volumes  are  given  by 


n,n+l 

f.  . 

1.3 


cartesian  ,  n,n+l  n,n+l,  ,  n,n+l  n,n+l, 

x-y  ■  V,  -  sra_i  >‘Vi  ‘Vi 


<  cylindrical  ..  n,n+L2  .  n,n+1.2,f  n,n+l  n,n+l, 

1  <Vi  ^  -  lVi  >  1  ‘Vi  -  Vi  1 


,  spherical  2-t  n  ,n+l  3  ,  n,n+l,3,  .  n,n+l  n,n+l, 

L  r-3  :T[(xi+i  >  -  (Xi-i  ’  1[c°s  yH  -  COS  yj+i  1 


We  have  then 


T 

„n+l  x  „n  n 

V.  .  o  .  .  =  ¥.  .  p .  . 

1,3  1,3  i,3  V3 


x..  n  „n+}  rtTn+i  . 

ot  (p  .  .  .A.  ,  .  6U.  ,  .) 

1-4,3  l+i ,  3  1-4,3 


n  »n+$  rT,n+i  .  n+J 

+  6t(p.  ,  ■.  A.  ,  .  <5U .  ,  .)  +  source  .  . 

i~t , 3  i~J » 3  1-1,3  i,3 


(332a) 
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and 


*n+1  p^.  =  ¥n  .  pn  . 

1.1  1.1  1.1  1.1 


...  n  *n+i  ,„n+4  % 

6t  (p  .  .  ,  A .  .  ,  £V .  .  , ) 

i,l+4  i.i+i  l.i+i 


+  6t(pn  .  ,  An+^  .  <5Vn+^  )  +  source  n+^ 

i.]-J  i,l-i  i.l-i  i.l 


. #N  )  ,  where 

y 


for  (i  =  )  and  (j  =  1,. 

n  .  n 
D  .  .  +  p  .  , 

0n  =  i±La 

Pi+4,j  2 


for  i  =  1,...,N  -  1  and  j  =  1,...,N  ,  while 

x  y 


„n  n 

n  _  ^ i , j  pi,j+l 
Pi,j+4  2 


for  i  =  l,...,Nx  and  j  =  1,...,N  -  1.  The  boundary  values  p 

are  obtained  from 

n  ^  n 

n  _  Pl'l  PL,i 
*  '  2 


4 ,  j 


Nx+4,i 


n  ,  n 
PNx,j  PR,j 


for  (j  =  1,...,N  )  where  L  and  R  denote  left  and  right  boundaries, 
respectively,  while 
n  n 

n  _  °i,l  Pi,B 

°i,4  ~  2 


n  n 

pi,N  +  Pi,T 


i»Vi 


for  (i  =  1,...,N  )  where  B  and  T  denote  bottom  and  top  boundaries, 


(332b) 


(33a) 


(333b) 

V>'i 


respectively.  The  mean  interface  areas  A*?*|  j  and  A^+^+^  are  given  by 


An+i  ,  n+i  _  n+i 
l+i t j  Yj+i  Yj~i 


(334a) 


. n+i  n+i 

A  ,  =  x  7 

i,j+i  i+i 


Xi4' 


(334b) 


where 


n  n+1 

n+i  yj+i  +  yi+i 


(335a) 


n  n+1 

,  x.  ,  +  x.  , 

n+i  r+i  r+i 


(335b) 


for  cartesian  x-y  coordinates;  by 


.n+i  „  n+i ,  n+i  n+i, 

‘  2"  WVi  -  yj-i 


(336a) 


An+i  . 
i/3+i 


u":*)2]. 


(336b) 


where 


n  n+1 

n+i  y1+i  +  yj+i 


(337a) 


r  (xn  }2+  n+1  2ni 

n+i  _  tXi+r  (  i+iJ 

Xi+i  “  2 


(337b) 


for  cylindrical  r-z  coordinates;  and  by 


'■  *■  net*** 


An+{  . 
i+i/3 


An+i  . 
i,]+l 


2*(x^j)2[cos  -  cos  y^t] 


.  n+J  2 
1lIxi+i) 


.  n+i, 2,  n+$ 

“W  ]  sin  yj+i. 


where 


r  1  r  .  _  n 


i  1  <  2  r  -L  r  **  'ITA  . 

yj+i  =  arc  cos  [-{cos  y.+i  +  cos  yj+i}] 


I-  ,  n  .3  ,  n+l  3-1 

[W  1 

2 


for  spherical  r-9  coordinates.  Finally, 


n+i  n+i  g 

OU.  I  .  =  u.  ,  .  -  u.  .  . 

i+i/3  i+i,]  i+i,3 


..  n+i  n+i  g 

6V.  .  ,  =  v.  .  ,  -  v:  .  , 

i/3+i  i/3+i  i/3+i 


The  grid  velocities  ug  ,  . ,  vg  .  ,  are  given  by 

i+i/3  i/J+i 


ug  .  . 

i+i/D 


vg  .  , 
i/3+i 


n+1  n 

x .  ,  .  -  x .  , 

i+j  i+i 


n+l  n 
*j+*  '  yj+* 


for  cartesian  x-y  coordinates;  by 


Ui+i  /  j 


,  n+l  2  .  n  .2 

-  (W 

2  xn+J  6t 
i+i 


(338a) 


(338b) 


(339a) 


(339b) 


(340a) 


(341a) 


(341b) 


(342a) 


•* 
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n+1  n 

vs  .  h ±Ll  ■j+i 

i ,  j+i  6t 


for  cylindrical  r-z  coordinates;  and  by 


,  n+1, 3  ,  n  , 3 

g  (Xi+j}  ~  (Xi+^ 

'i+i'j  “  3(x£j)2  St 


vg  . 


i,j+i  5t 


n+i  n  n+1 

x.  cos  y  .  , ,  -  cos  y . , , 
1  3+* _ lit! 


sm  y 


n+i 

j+i 


for  spherical  r-0  coordinates,  where 

.  n+i  2  ,  n+i  n+i  n+i  2 

,,  ->  (x.  ,)  +  (x  ,)  (x.  .)  +  (x.  ,) 

n+i  _  2_  i+j _ i+i  i-j _ i-j 

xi  3  n+i  n+i 

x  ,  +  x  , 

i+i  i-i 

Equations  (334)-(344)  are  valid  for  i  =  0,  1,...,!*  and  j  =  0,  1,...,J 

Tx  Ty 

Equations  (332)  yield  p.  .  and  p.  .,  which  are  used  later  to  evaluate 

x ,  3  x ,  j 

the  antidiffusion  fluxes.  The  transported  and  diffused  densities  are 
obtained  from 


.  n+1  TD  „n+l  T  „n+l  n  n 

¥.  .  p .  .=¥  .  p .  .  +  v .  ,  .  V.  ,  . (p .  ,  .  -  p .  .) 

1,3  X,]  1,3  x,3  x+i,3  l+i,3  1+1,3  1,3 


TIn+l  ,  n  n  ,  ,  ,  n+1  ,  n 

-  V.  -  .  ,  . (p .  .  -  p.  ,  .)  +  V.  .  ,  ¥.  .  ,(p.  .  . 

x-i,3  x-i,3  x,3  x-1,3  i,3+i  x,3+i  i,3+l 


n  ,  „n+ 1  ,  n  n  , 

-p.  . )  -  v .  .  ,  ¥.  .  ,  (p .  .  -  p  . 

i,3  i , 3~i  x,3-i  x,D  i,3“l 


for  i  =  1 , . . . ,N  and  j  =  1,...,N  ,  where 
x  J  y 

_  1  12 

Ji+i , j  6  +  3  ei+i , j 


(342b) 


(343a) 


(343b) 


(344) 


l  • 

y 

then 


(346) 


(347a) 


and 
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_  1  12 

./ j+i  6  +  3  £i, j  +  i' 


while 


v"+1  =  I(^+1  +  vn+1  ) 

vi+i,j  2^i,j  i+l,  j' 


for  i  =  1,...,N  -  1  and  j  =  .  Similarly, 


V.  «  i(Vn+1  +  ¥n+1  ) 

i,3+i  2  1,3  i,3+1 


for  i  =  1,...,N  and  j  =  1,...,N  -  1.  At  the  boundaries, 

x  y 


n+l  „n+l  ,  n+1  _  n+1 

^  N  +},j  ¥N  ,j 


for  j  =  1, .. . ,N  ,  while 


V")  -  and 

i/i  x,l  x / Ny ^*2  y 


for  i  =  1,...,N  .  The  dimensionless  velocities  ei+^j, 
obtained  from 


are 


5U .  AntJ  .6t  .  , 

=  i+},j  i+},3  ,  1  + _ L 

£i+i,j  2  vn+l  ^n+1 


-) 


i ,  j  i+l, j 


for  i  =  0,...,N  and  j  =  1,...,N  using  (349a)  and 
x  y 


£i, j+} 


,„n+i  _n+} 

oV .  *,1  A .  ,  i  vt 

1,3+i  i,3+} 


¥n+1  ¥n+1  1 
i,3  1,3+1 


for  i  =  1,...,N  and  j  =  0,...,N  using  (349b).  The  antidiffusion 
x  y 

are  then  evaluated  according  to 


-r 


(347b) 


(348a) 


(348b) 


(349a) 


(349b) 


(350a) 


(350b) 


fluxes 
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In  the  above  p^*,  are  fc^e  uPPer  an<^  lc>wei  bounds,  respectively, 

on  pn+^,  chosen  so  as  to  guarantee  no  ripples  formation  at  grid  point  (i,j). 
i  / 1 

Finally,  since  each  flux  leaves  a  cell  to  enter  an  adjacent  one, 

3.  The  fluxes  correction  factors  are  defined  as 


Ci+4,j 


^min  (R+  ,  .,  R.  , )  if  F .  ,  .  >  0 

1+1,3  1,3  1+4,3  ~ 


Lmin  (R,  ,  .,  R.  .)  if  F .  .  .  <  0 

1+1,3  i,3  1+4,3 


(356a) 


and 


Ci, j+4 


min  (R.  .  , ,  R.  . )  if  F .  .  ,  >  0 

1,3+1  1,3  1,3+4  - 


*min  (R.  .  ,  ,  R .  . )  if  F .  ,  .  <  0. 

1,3+1  1,3  1,3+4 


(356b) 


The  corrected  fluxes  are  given  by 


pc  =  C  F 

i+4 , j  i+4 , j  i+4,j' 


(357a) 


Fi , j+4  Ci,j+4  Fi , j  +  4  ' 


(357b) 


,  _  max  ,  min  .  . 

4.  For  p and  p ,  two  choices  are  presented.  A  conservative 


choice  would  be 


max 

D  .  . 

1.3 

min 
p  .  . 

1.3 


TD 

TD 

TD 

TD 

TD 

i-1,  j' 

Pi,j-1' 

Pi,3' 

Pi+l,j'  Pi , j+l 

TD 

TD 

TD 

TD 

TD 

i-1, j' 

P  ■  -i  t 

1,3-1 

Pi,j' 

Pi+l,j'  Pi , j+l 

)  ; 

)  . 


(358a) 


(358b) 


A  more  tolerant  choice  that  gets  rid  of  the  problems  of  "clipping"  and 
"terracing"  partially  is 


max  ,  a  a  a  a  a 

p.  .  =  max  (P .  .  .,  P.  .  ,,  P.  .  ,  p...  .,  P. 

1,3  1-1,3  i,3“l  1,3  1+1,3  1,3+1 


(359a) 
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where 


^a 

"i,  j 


.  TD  n  . 
max  (p  .  .  ,  p  .  . )  , 

1,3 


and 


min 


min 


.  / 

-3 


,b  b 

"i,j-l'  *"  i ,  j  ' 


l+l ,  j 


where 


(359b) 


b 

o .  . 

i-3 


,  TD 

mm  (p  .  . 

i-3 


n  , 

P  .  . )  . 

i-3 


ANTIDIFFUSION  AND  HALF-STEP  UPDATING 

The  corrected  antidiffusion  fluxes  are  added 


n+1  n+1  _,n+l  TD  ,_c  _c  , 

V.  .  c  .  .  =  V.  .  p .  F.  .  .  -  F.  ,  •)  - 

1-3  1-3  1-3  1-3  1+4-3  1-4,3 


thus  giving  the  new  density  o  . 


n+1 

i-j 


n+4 


Finally,  it  remains  to  specify  u  .  . 

i+i  -  3 
n+  4 

and  the  source  terms  denoted  bv  "source 

i-3 

interfaces  are  obtained  from 


(F 


c 

i .  j  +  4 


(360) 


,  n+4 

and  v .  .  . 

i.3+4 

First,  the 


in  equations  (340) 
velocities  at  the 


n+4  n+4 

,  u.  .  +  u.  , 

n+4  =  _ i+iil 

i+4-j  2 


(361a) 


for  i  =  1,...,N  -  1  and  j 

x 

n+4 

,  u .  .  +  u. 

n+l  _  i-3  u 

4-j  ~  2 


=  1 


n+4 

u„,  .  +  u„ 

.  N  ,  3  R 

n+4  _  x  _ 

+4,j  2 

x  J 


,N  while 

y 


and 


n+i  n+i 

v .  +  v  .  , 

i,l  i,i+l 


i,  j  +  i 


(361b) 


for  i  =  and  j  =  1,., 

n+i 

n +  i  Vi,l  +  VS 
Vi,i  ~  2 

n+i 

,  Vi,N  VT 

vn+*  =  — y _ 

i,Ny+i  2 


,N  -  1  while 

y 


As  for  the  source  terms  they  were  defined  earlier,  Eqs.  (321)  through  (329) 
.  .  n+i  n+i  ,  n+i 

Next,  to  get  u.  . ,v.  .,  and  source  . ,  we  advance  our  system  of  conserva- 

1,3  1,3  i,3 

.  .  ......  n  n  n 

tion  equations  i  time  step  using  u.  .,  v.  .,  source  .  .,  tnen 

i , j  i,j  i , 3 

n+i  i  n  n  ..  n+1  j  n  n  ot 

u .  t+t  +  ot  =  u .  .  t+t  +  — 

i,3  1  i,3  1  2 


n+i  I  *n  n  .  r*.  n  +  1  I  _n  n  3t 
V .  .  t+t  +  6t  a  V .  .  t+t  +  -r~ 
1,3  i,3  2 


n+i  |  ,  n  n  .  n+1  i  .n  n  ct 

source  .  .  t+t  +  ot  =  source  .  .  t+t  +  — 
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XVII.  TWO-DIMENSIONAL  TIME  SPLITTING 
VERSUS  FUL  v  TWO-DIMENSIONAL  ALGORITHMS 


Going  back  to  Eq.  (202) , 

A(3x,3y)  =  e  lUxBx  +  £yBy>  =  e  1£x6x  e  1'yBy  =  A(3x)A(6^)  (362) 

If  P .  .  =  e  ,  where  x  =  (i6x,j5y),  the  analvtic  solution  of  t—  +  u-V;  =  C, 

i,2  '  at 

according  to  Eq.  (362),  yields 


■  AivM6x>oLi 


where  u  =  (u,v)  is  constant  and  the  two  operators  A(8x)  and  A(8^)  are 
commutable.  Noticing  that 


x  . ,a  ,  n 
°i.j  "  A(8x)pi,j 


(364a) 


is  the  analytic  solution  of  4^-  +  u  =  0,  whereas 

dt  dX 


n+l  ,  ,  x 

’i.j  * 


(364b) 


"ho  x 

is  the  analytic  solution  of  —  +  v  ^  =  0  for  an  initial  density  p .  . , 

3t  ay  i,J 

equation  (363)  invokes  time  splitting  as  an  exact  solution  for  the 

linear  PDE.  If  we  derive  a  numerical  scheme  by  expanding  A(3X»8  )  in  terms 

of  sin  6  ,  cos  3  ,  sin  8  and  cos  8  such  that  both  A (8  ,8  )  and  its 
x  x  y  y  x  y 

expansion  agree  up  to  a  prescribed  order  of  S  and  6  ,  we  obviously  end 

x  y 

with  a  time  splitting  scheme,  in  which  each  of  the  x  and  y  operators 

agrees  with  A (8  )  and  A (8  )  up  to  a  prescribed  order  of  8  and  6  , 
x  y  x  y 

respectively . 

Alternatively,  if  a  1-D  scheme  is  n  -  order  in  phase  error 
and  m  -  order  in  diffusion  error,  namely 


*'  -  "  ,  -  -VC 


■'f  *  jt 
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!a|  =  i  +  o(g'n) 


(365) 


and 


e  =  e  +  o(sn+1) 

exact 


(365 b) 


where  |A|  and  3  are  the  amplitude  and  angle  of  the  scheme  transfer  function  A, 

i.e.,  A  =  [A |  ei9,  using  a  time-splitted  version  of  the  one-dimensional  scheme 

to  solve  a  two-dimensional,  x-y  problem,  gives  j  A  [  e1"'  =  A  S  A^Ay 

=  ( |  A  |e1,'x)(jA  j  e  y )  .  Thus, 
x  y 

|A|  =  |A  |  •  |A  |  =  (1  +  0(Sm))(l  +  0(Bm))  =  1  +  0(£m)  +  0(Sm)  (366a) 

x  ’y  x  y  x  v 

and 

3  =  5  +  3  =  [5  +  o(3n+I)l  +  [3  _  +  0(3n+1)l 

x  y  x  exact  y  exact  y 


«.  ,_n+l.  ,_n+l. 

'exact  +  0(Bx  >  +  0(3y  ) 


(366b) 


showing  the  two-dimensional  scheme  to  be  of  the  same  order  as  the  one¬ 
dimensional  one.  Moreover,  the  errors  in  both  j  A  j  and  3  are  free  from 

nl  n2 

mixed  freauencies,  such  as  n(6  3  )  where  n.  +  n_  =  m  or  n  +  1. 

x  y  12 

Although  time-splitting  appears  to  be  the  perfect  solution, 
physically  unacceptable  results  are  produced  when  dealing  with  incom¬ 
pressible  or  nearly  incompressible  flow  fields,  or  when  a  differential 
identity,  such  as  divergence  free  property  or  irrotationality ,  is  to  be 
strictly  enforced.  Moreover,  because  the  antidiffusion  fluxes  are 
corrected  in  each  direction  independently  of  the  other,  unnecessary 
"clipping"  occurs.  Namely,  the  flux  corrector  may  cancel  a  flux  that 
would  produce  a  ripple  in  one  direction,  which  actually  is  safe  in  two- 
dimensions  due  to  the  growth  or  decay  of  the  adjacent  cells  in  the  other 
direction . 


vr  ry 


Going  back  to  the  problems  arising  in  incompressible  flows,  let's 
consider  a  case  where  U  =  U(x,y),  independent  of  time,  satisfying  7-U  =  0 
For  simplicity  assume  u  =  uq  +  Cx  and  v  =  vq  -  Cy.  Figure  41  illustrates 

the  velocities  at  the  interfaces  of  a  cell,  when  5x=6y=l,C=0.1. 

* 


0.1 


0.1 

\ 

if 

• 

3 

£  7 

=1 

„  5x=l 

l _ 

0.2 

- > 


0.2 


Fjg.  41 


Using  a  simple-transport  scheme  with  time  splitting 


V.  . 
1,3 


¥.  .  p1  .  =  V.  .  px  . 

X,:  Kx,3  x , 3  px,3 


,  o  o  ,  .  „ 

(u,  .  p .  ,  .  -  u.  .  .  ; .  ,  .) :y6t 

x+i,3  x+i , 3  x-i , 3  i~i,3 


(v.  .  ,  pX  .  .  -  v.  .  ,  pX  .  ,)£xft 

i,3+i  i/D+i  i,3-i  x»3“t 


(367a) 


(367b) 


where  0,1  stands  for  t  =  0,  it,  respectively.  Assuming  a  uniform  initial 

density  o°  =  1,  and  ot  =  1,  Eq.  (367a)  gives  (1) •  (pX  ,)  =  (l)-(l)  -  ((0.2) •  (1) 

1 » 3 

-  (0. 1) • (1) ) • (1) • (1)  yielding  p.  .  =  0.9.  Since  u  =  u(x),  v  =  v(y),  pX  .  =  0. 

1 ,  3  1 »  3 

for  all  j's  and  since  u,v  are  linear,  it  is  also  true  for  all  i’s.  Then,  from 

Eq.  (367b),  we  obtain  (l)-(p1  .)  =  (1) -  (0.9)  -  (  (0.1) -  (0.9)  -  (0. 2) •  (0. 9) ) • 

1 ' 3 

(1) •  (1)  yielding  p.  .  =  0.99.  After  n  time  steps,  it  is  obvious  that 

i,3 

oa  .  =  (0.99)n  for  all  i  and  j.  Generally,  p°  .  =  p°  .  (1-C“)n.  In  other 

1,3  i<3  x,3 

words,  the  density  keeps  on  uniformly  distributed  but  decreases  with  time 


continuously . 
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An  equivalent  fully  two-dimensional  scheme  would  be 


1  o  o  o  ,  .  . 

V.  .  p.  .  =  V.  .d.  .  -  (u .  .  .  o .  ,  .  -  u .  .  .  p .  ,  .  jyot 

1,3  1,3  i/l  1/3  1+4,3  1+4,3  i-J/3  i-i/ 3 


,  o  o  .  c  P 

(v .  .  ^ ,  P .  .  -  v .  .  ,p  .  .  .  6x5 1 

1,3+4  i+4  1/3-4  1,3-4 


(368) 


which  gives  p.  .  =  1,  i.e.  conserves  the  mass, 
i  /  3 

The  discrepancy  obviously  lies  in  the  assumption  of  U  =  const 
while  p  is  varying  when  deriving  Eq.  (362) .  In  terms  of  transfer  func¬ 
tions,  the  scheme  of  Eqs.  (367)  is  written  as  A  =  (1  -  e  t  )  (1  -  e  t  ) 

x  x  y  y 

whereas  that  of  Eq.  (368)  takes  the  form  A  =  1  -  s  t  -et. 

x  x  y  y 

The  difference  is  obviously  in  the  term  "sett  "  which  as  will  be 

x  y  x  y 

shown  later  is  essential  for  high  order  diffusion.  In  the  next  section,  we 
try  to  cast  a  time-splitted  scheme  into  a  fully  two-dimensional  form.  A 
detailed  explanation  of  the  problems  involved  is  given. 


FULLY  TWO-DIMENSIONAL  VERSIONS  OF 
TIME-SPLITTED  SCHEMES 

Going  back  to  the  fourth  order  phase  and  diffusion  scheme 


A  =  (1  -  et)  (1  -  yd) 

1  c2 

where  v  =  —  +  —  and  y 
6  3 

of  Eq.  (370) 


+  vd 


6 


(370) 

The  two-dimensional,  splitted  version 


((1  -  extx)U  - 


V  d  ) 
x  x 


W*[<i  - 


s  t  )  (1  - 

y  y 


y  d  ) 

y  y 


v  d  ] 

y  y 


(371a) 


or 
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A  =  (1  -  e  t  +  v  d  )  (1  -  e  t  +  v  d  )  -  y  d  (1  -  et  )  (1  -  c  t  +  v  d  ) 

XX  XX  y  y  y  y  XX  y  y  y  y 


-yd  (1  -  e  t  )  (1  -  e  t  +  v  d  )  -  y  y  d  d  (l-et)(l-st> 
yy  yy  xx  xx  xyxy  xx  y  y 


(371b) 


can  be  written  as 


„  etvd  etvd 

A™  -  1  -  e  t  (1  -  -P  *  -P)  -  £  t  (X  -  -f-i  ♦  -P) 
xx  2  2  yy  2  2 

et  vd  et  vd 

,  y  y  y  y,  ,  xx  xx, 

+  v  d  (1 - X-X-  +  )  +  v  d  (1 - - —  +  — — ) 

xx  2  2  yy  2  2 


(372a) 


which  is  >0  for  [e  |,!e  I  <  i,  therefore  ensuring  positivity  of  o  if 


x '  1  y  1  — 


0  >0.  Then, 


TD  *  1 

A  =  A  -yd  (1  -  e  t  )  [1  -  e  t  +  v  d  --gyd(l-et)] 
x  xv  x  x  L  yy  yy  2  *y  y  yy 


u  d  (1  -  e  t  )  [1  -  £  t  +  vd  -  ^  |i  d  (1  -  e  t  ) 

yy  yy  xx  x  x  2  x  x  xx 


(372b) 


where  the  asterisks  denote  the  operators  which  fluxes  are  to  be  corrected. 
This  will  allow  us  to  correct  the  x  and  y  antidiffusion  fluxes  simul¬ 
taneously,  thus  avoiding  unnecessary  clipping.  But  that  does  not  solve  the 
problems  associated  with  divergence  free  flow  fields,  for  example,  because 

of  the  term  "s  t  t  t  . "  Moreover,  we  notice  that  the  form  of  Ea.  (372)  is 
xyxy 

in  no  way  unique. 

Although  Eqs.  (366)  show  in  a  clear  simple  way  that  A  =  A^A^ 

fourth  order  in  phase  and  diffusion,  let  us  analyze  it  using  Eqs.  (246)  to 

(243),  Eqs.  (260)-(262)  with  Eqs.  (241)-(244).  The  purpose  is  to 

determine  which  terms  are  responsible  for  the  fourth  order  diffusion, 

fourth  order  phase  error,  positivity,  stability,  and  so  on.  We  notice 
t  t  t  t 

that  A  X  X  =  0,  A  ^  y  =  0,  making  Eqs.  (241) -(244)  valid. 


106 


Differentiating  Eq.  (371a) ,  we  get 


x 

A  y  =  -  £  (1  -  y  d  )  A 

X  X  X  V 

y  y  y  * 

d 

X 

A  y  =  [ (v  ■  u  )  +  e  y  t  ]A 
x  x  xx  x  y 

y  y  y  y  y  * 

t  d 

X  X 

A  y  Y  =  £  y  A 
x  x  y 

y  y  x 


A  y  =  ■  y1  ■  w*vl(vx  “  V +  YxV 


1 1 

A  x  y  =  e  e  (1  -  yd)  (1  -  yd) 
x  y  xx  y  y 


1  *  y  *  1  -  “»>  +  VxV 1  (”y  ■  V  *  eygyty! 


(378a) 


(378b) 


At  3x  =  6y  =  0,  tx  =  ty  =  0,  and  dx  =  dy  =  0,  thus  Ax  =  Ay  *  1 ,  yielding 


°  ? 


y 

A  =  v  -  y„ 

o  XX 

y  y 


y  y 

A  x  *  =  £  y 

O  XX 


aq  =-yvv 


x  y 

A  *  *  £  € 


d  d 

A  x  y  =  (v  -  y  )  (v„  -  y„) 


(383a) 


(383b) 


x  x  y  y 


*  -  -  •  »  * 
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Substituting  into  Eqs.  (246)  and  (248)  ,  we  get 


x 

(log  A) ^  =  -  ie 
o  x 

y 


(384) 


(log  A) 


XXX 

yyy 


it  [1  -  6  (v  -  u  )]  -  i[6e  u 

X  XX  XX 

y  y  y  y  y 


2£J1 

X 


1 

=  6ie  (-  +  -f  -  V  ) 
x  6  3  x 

y  y 


(385) 


showing  A^A^  to  he  fourth  order  in  phase  error,  but  more  importantly,  that 
the  cross  terms  of  Eqs.  (382)  and  (393a) ,  which  do  not  appear  in  one- 
dimension,  are  essential  to  the  fourth  order  phase.  More  specif ica 1 ly , 
these  cross  terms  reduce  the  dependence  of  phase  error  on  v  and  u  to 
dependence  on  v  only,  leaving  y  free  to  be  adjusted  for  a  high  order 

diffusion. 

d  d 
x  v 

A  in  Ea.  (383b)  is  not  used  in  either  (384)  or  (385)  and 

o 

therefore  can  take  any  value  without  affecting  the  phase  error. 

Now  we  can  construct  the  simplest  fourth  order  phase  error  scheme. 
Such  a  scheme  has  to  satisfy  Eqs.  (379)  to  (382)  plus  Eq.  (383a)  giving, 


A=(l-et)(l-et)  +  (v  -y)d  +(v  -  t  )dy  *  c  t  u  d 

xx  yy  x  xx  y  y  xxxx 


+  e  t  u  d  -et(v  -  y  )d  -  e  t  (v  -  u  )d 
y  yy  y  xxy  yy  yyx  xx 


(386) 


where  the  integration  constant  was  selected  as  unity  to  satisfy  consistency, 
i.e.,  A(3x,6y  =0)  =1.  Eq.  (386)  can  be  written  as 

A  =  (l-et)(l-et)  +  vd  (1-et)  +  v  d  (1-et) 
xx  yy  xx  yy  yy  xx 


-yd  (1-et  -  et)  -  yd  (1-et  -  et) 

xx  xx  yy  yy  xx  y  y 


(3S7) 


d  d 

Since  A  X  ^ 

o 
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does  not  affect  the  phase  error  order,  we  can  assign  a  value 

for  it  that  would  ensure  positivity.  We  add  to  the  terms  of  Eq.  (387) 

"v  v  d  d  yielding 
x  y  x  y 

A  -  (1  -  e  t  +  v  d  )  (1  -  e  t  +  vd)-ud(l-et  -  e  t  ) 
xx  xx  yy  yy  xx  xx  yy 


-  p  d  (1  -  e  t  -  et) 
y  y  XX  y  y 


Now,  substituting  Eqs.  (379)  to  (383)  into  Eqs.  (242)  and  (243)  ,  we  get 
x 

Ay  =  -  ie  C 

o  x 

y 


2 (v  -  u  ) 

X  X 

y  y 


(390a) 


xy 

A  =  -  z  t 


(390b) 


which  when  substituted  into  Eqs.  (261) ,  yield 


2  xx 

(IaI  )YY  =  2[-2(v  -  u  )  +  e  ]  -  0 

1  1  O  XX  X 

y  y  y 


(391a) 


‘W2^-21-  Vy  +  exJ  =° 


(391b) 


showing  A  A  to  be  fourth  order  in  diffusion  error, 
x  y 

t  t 
x  y 

Notice  that  A  =  s  £  ,  is  essential  for  fourth  order  diffusion 

o  x  y 

(already  satisfied  by  the  scheme  of  Eq.  (388)). 

The  simplest  fourth  order  (phase  and  diffusion  error)  positive 
scheme  is  therefore  that  of  Eq.  (388).  It  is,  however,  unstable.  For 
instance , 


(  -  "  ♦  - -T  ■■ 
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A 


R 


1  +(v 

x 


)d 
x 


U  )  d 

y  y 


vvdd  +  e  e  t  t 
xyxv  xyxy 


(392a) 


while 


A  =  -  £  t  [1  -  u  d  +  (v  -  u  )  d  ]  -  e  t  [  1  -  v  d  +  (v  -  |j  )d  ]  (392b) 

I  xx  xx  y  yy  yy  y  y  x  x  x 


At  e  =  e  =  i  and  6  =  g  =  tt/2,  d  =  d  =  -2  and  t  =  t 

xy  xy  xy  ,  xy 


„  .  1111, 

AR  1-4_4  +  4~4 


--  ill  *i-|l  -  ill  .i-il  ■  -  1 


i,  yielding 

(393a) 

(393b) 


Since  we  know  that  A  of  Eq.  (372)  is  stable,  let's  try  to  approach 
it  in  steps.  First,  we  try 


A  =  (1  -  c  t  +  v  d  )  (1  -  s  t  +  vd)-vid(l-et)(l-£t) 
xx  xx  yy  yy  xx  xx  yy 


-ud(l-et)(l-et) 
yy  y  y  x  x 


(394) 


thus  adding  "-(u  d  +  u  d  ) (c  £  t  t  ) "  to  the  real  part,  becoming  then 
xxyyxyxy 


A  =  1/2  -  1/8  =  3/8 
R 


(395) 


still  unstable.  Next,  we  try 


A  =  (1  -  e  t  +  v  d  )  (1  -  £  t  +  v  d  )  -  u  d  (1  -  £  t  )  (1  -  £  t  +  v  d  ) 

xx  xx  yy  yy  xx  xx  yy  yy 


pd(l-et)(l-et  +  vd) 
yy  v  y  xx  xx 


(396) 


This  will  add  ”-(y  v  +  u  v  )d  d  "  to  the  real  part  and  "(e  t  y  v  ) 
x  y  y  x  x  y  x  x  x  y 

+  eytylJyVx^dxd^"  to  the  imaginary  one.  We  get  then 


i 
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(397a) 


A  =1-A=± 
R  8  4  8 


AI  = 


(397b) 


-j  ^  1  2  50 

whence  J  A ( *"  =  (— )  +  (— )  =  —  c  1,  showing  Eq.  (396)  to  be  too  stable  at 

3  =6  =  m  .  Moreover,  it  is  not  phoenical;  Ablate  =  e  =0.  These 
x  y  2  x  y 

two  effects  can  be  avoided  by  picking 


A  =  (1  -  e  t  +  v  d  )  (1  -  £  t  +  vd)  -  u  d  (1  -  e  t  )  (1  -  e  t  +  —  v  d  ) 
xx  xx  yy  yy  xx  x  x  yy  2yy 


ud(l-st)(l-et  +  —  v  d  ) 
yy  yy  xx  2xx 


(398) 


yielding 


n  -  3  1  1  _  1 

AR  8  “  2  X  4  4 


_  .11  15 

A  =-l+~x— ■=  -  — 
I  2  8  16 


(399a) 

(399b) 


2  ^5  2  ^  2  ^41 

whence  ( A  j  =  (— )  +  (— )  =  <  1,  closer  to  1,  therefore  promising  a 

smaller  net  diffusion  and  phoenical  since  A  =  1  at  e  =£  =0. 

x  y 

Noticing  that  the  added  terms  to  Eq.  (388)  are  triple  operators 

1 1 1 

(ttd  ,  ttd  ,  tdd  ,  t  d  d  )  they  have  no  effect  on  (log  A) 
xyxxyy  xxyyxy  J  o 

Eq.  (398)  is  still  fourth  order  in  phase  error.  Furthermore,  upon  expanding 
j  A  j  2 ,  we  get 


!  A 1  2  =  1  +  ~—[  e2  (1  -  e2)S4  +  [£2(  1  -  £2) 

11  12  x  xx  x  y 


2  2  2  2  2  2  4 

+  £  (1  -  £  )]8  6  +  £  (1  -  £  )6  }  +  ... 

y  x  x  y  y  y  y 


(400) 


showing  diffusion  error  to  be  of  fourth  order.  The  scheme  is,  however, 


slightly  unstable  near  8  =  8  =0,  since  the  fourth  order  coefficient  is 

x  y 


Ill 


2  2 

positive.  We  notice  also  the  presence  of  a  term  "3^6^"  in  Eq.  (400)  (also 
in  the  phase  error  expansion) ,  which  does  not  show  in  the  expansion  of  Eq. 
(372)  (according  to  Eq.  (366) ,  making  the  scheme  of  Eq.  (398)  slightly 
inferior  to  that  of  Eq.  (372) . 

Upon  comparing  Eq.  (398)  to  (372) ,  it  is  obvious  that  Eq.  (372) 

cannot  be  much  simplified;  at  least  without  sacrificing  stability  or 

phoenicity.  Whichever  we  use,  the  n°  of  operations  involved  in  evaluating 
n+1 

o  is  much  larger  than  that  in  the  fully  two-dimensional  scheme  of  Eq. 
(226e) .  Moreover,  the  n°  of  two-dimensional  arrays  required  to  store  the 
intermediate  values  is  enormous. 

Since  the  only  advantage  of  Ecs.  (372)  ,  the  fully  two-dimensional 
version  of  the  time-splitted  scheme  of  Eq.  (371b)  is  the  reduced  clipping 
associated  with  the  flux  limiter,  we  conclude  that  time  splitting  is  the 
sensible  answer.  We  abandon,  therefore,  trials  to  cast  the  time-splitted 
scheme  in  fully  two-dimensional  versions. 


»  -  -t  *- 
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X VIII .  IMPROVING  DIFFUSION  ERROR  OF  THE 


FULLY  TWO-DIMENSIONAL  SCHEME 


Now  that  we  have  classified  the  terms  responsible  for  the  fourth 

order  phase,  diffusion,  etc.,  in  the  time-splitted  scheme,  let's  go  back 

to  the  fully  two-dimensional  scheme  and  study  the  terms  preventing  us 

from  reaching  a  fourth  order  diffusion  error.  As  explained  earlier,  the 

term  "sett  "  is  essential  to  reduce  the  dependence  of  the  phase  error 
x  y  x  y 

to  one  on  v  alone,  thus  leaving  y  free  to  be  adjusted  for  a  high  order 

diffusion.  A  closer  look  reveals,  however,  that  the  above  conclusion  is 

an  indirect  one.  The  direct  conclusion  is  that  "c  e  t  t  "  is  needed  to 

x  y  x  y 

cancel  3  "  resulting  from  squaring  the  imaginary  part.  Specifically, 

any  scheme  has  to  incorporate  the  combination  (e  t  +  e  t  )  leading  to 

xx  y  y 

its^  sin  3x  +  sin  3^ )  which  is  approximated  by  i/s^B^  +  £^3^).  To  car>cel 

it,  a  term  including  sin  3  sin  3  is  needed.  Besides  t  t  ,  the  above 

x  y  xy 

term  can  also  result  from  cos  (3  :5  ),  i.e.  diagonal  diffusion 

x  y 

(pn  ...  -  2cn  .  +  cn  ,  .  ,).  Admitting  diaaonal  terms  is  outside  the 

l+l, ji.'  1,3  l-l, 1+1 

scope  of  this  article  and  is  left  out  to  an  upcoming  one.  However,  we 

empnasize  that  there  is  a  stability  problem  caused  by  the  imaginary  part 

iA_  =-[st  (l  -  u  d  )  +  c  t  < 1  —  u  d  ) ]  which  amplitude  is  already 
I  xx  x  x  v  v  v  y 

1  TT 

larger  than  unity  for  c  =  £  =  — ,  8  =3  =  m  unless  y  =  y  =0 

xy2xy2  xy 

there,  in  which  case  we  have  a  large  residual  diffusion.  Adding  just  a 
diagonal  diffusion  can't  help,  since  it  only  adds  to  the  real  part. 
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APPENDIX 


SUBROUTINE  FCT2D (RHO, KO . KN, KR, SNKRNZ , 

&  LBC . RHQLBC ,  RBC ,  RHORBC , BBC , RHOBBC >  TBC , RHOTBC ) 

C 

C 

C  ORIGINATOR  :  RAAFAT  H.  GUIRGUIS 

C 

C  DESCRIPTION  : 

C  - 

C 

C  A  FULLY  2-D  ROUTINE  THAT  SOLVES  GENERALIZED  CONTINUITY  EQUATIONS 

C  OF  THE  FORM 

C 

C  D  RHO  /  D  T  =  -  DIV<  RHO  *  V  )  -  SOURCES 

C 

C  WHERE  RHO  IS  THE  GENERALIZED  DENSITY,  AND  V  IS  THE  FLUID  VELOCITY. 

C  FOR  SECOND  ORDER  ACCURACY,  IT  IS  ADVISABLE  TO  ADVANCE  HALF  A 

C  TIME  STEP  USING  THE  VELOCITY  AND  SOURCE  TERMS  AT  THE  BEGINNING 

C  OF  THE  TIME  STEP,  THEN  ADVANCE  A  WHOLE  TIME  STEP  USING  THE 

C  HALF-POINT  VELOCITY  AND  SOURCE  TERMS.  USING  THE  HALF  POINT 

C  DENSITY  IS  NOT  RECOMMENDED.  IT  IS,  HOWEVER,  INCLUDED  AS  AN 

C  OPTION,  BY  ALTERNATING  (KO, KN)  BETWEEN  (1,2)  FOR  THE  HALF  TIME 

C  STEP,  AND  (2,1)  FDR  THE  WHOLE  TIME  STEP. 

C  THE  OLD,  (KO) ,  AND  NEW,  (KN) ,  DENSITIES  (  AT  THE  BEGINNING  AND  END 

C  OF  THE  TIME  STEP,  RESPECTIVELY  )  ARE  STORED  IN  A  2-LEVEL  2-D  ARRAY 

C  (  3-D  ARRAY  ).  THE  MASS  AND  DIFFUSION  FLUXES  ARE  EVALUATED  USING 

C  KO  DENSITY,  WHEREAS  KN  DENSITY  DETERMINES  THE  ANT I -DIFFUSION 

C  FLUXES.  IT  IS  ADVISABLE  TO  SET  KO  =  1,  KN  =  2,  UNLESS  THE  HALF 

C  POINT  DENSITY  IS  TO  BE  USED  DURING  THE  WHOLE  TIME  STEP.  THEN 

C  KO  =  2,  KN  =  1  FOR  THE  WHOLE  TIME  STEP. 

C  KR  DETERMINES  THE  LOCATION  OF  THE  RESULTING  DENSITY.  IT  IS 

C  ADVISABLE  TO  SET  KR  =  2  ,  1  ,  FOR  THE  HALF  AND  WHOLE  TIME  STEPS, 

C  RESPECTIVELY.  THIS  CHOICE  ELIMINATES  THE  NEED  TO  COPY  THE  NEW 

C  DENSITY  ON  THE  OLD  ARRAY,  IN  PREPARATION  FOR  A  NEW  TIME  STEP. 

C  SNKRNZ  IS  A  LOGICAL  VARIABLE  WHICH,  WHEN  SET  TO  .TRUE.,  TELLS  THE 

C  ROUTINE  TO  USE  THE  CORRECTION  FACTORS  OF  THE  LAST  SNKRNZ  =  .FALSE. 

C  CALL,  TO  LIMIT  THE  ANT.I -DIFFUSION  FLUXES.  IF  SET  TO  .FALSE.,  THE 

C  CORRECTION  FACTORS  ARE  EVALUATED  FROM  THE  CURRENT  VARIABLES  AND 

C  USED  IN  THE  LIMITING  PROCESS. 

C  LBC,  RHOLBC ,  RBC,  RHORBC,  BBC,  RHOBBC,  TBC,  RHOTBC ,  ARE  DEFINED 

C  BELOW. 

C 

C  (1)A  PARTICULAR  GEOMETRY  IS  SELECTED  BY  A  CALL  TO  ENTRY 

C  SETGOM  : 

C  CALL  SETGOM (  4HCART ,  1HX,  1HY ,  NX,  NY  ) 

C  OR  . . , 1 HX ,  1HZ , . . .  OR  ...,  1HZ ,  1 HY ,..., ASSUMES  CARTESIAN 

C  COORDINATES.  ORDER  OF  THE  2  COORDINATES  IS  IMMATERIAL  ONLY 


c 

FOR  THIS 

CASE. 

c 

CALL 

SETGOM ( 

3HCYL , 

1  HR , 

1  HZ ,  NX 

,  NY  ) 

c 

CALL 

SETGOM ( 

3HCYL , 

1HR , 

3HFYE , 

NX,  NY  ) 

c 

CALL 

SETGOM ( 

3HCYL , 

1  HZ  , 

3HFYE , 

NX,  NY  ) 

C  FOR  THE  3  TYPICAL  CYLINDRICAL  COORDINATES. 

C  CALL  SETGOM (  3HSPH ,  1HR,  4HCETA,  NX,  NY  ) 


A1 


k-r  r— 


kZ 


.  ». 


Ur 


jt 


r 


i 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

n 

u 

c 

c 

c 

c 

c 

c 

c 

c 

c 


CALL  SETGOM (  3HSPH,  1HR,  3HFYE ,  NX,  NY  ) 

CALL  SETGOM (  3HSPH ,  4HCETA,  3HFYE ,  NX,  NY  ) 

FOR  THE  3  TYPICAL  SPHERICAL  COORDINATES. 

NX,  NY,  ARE  THE  NUMBERS  OF  CELLS  CENTERS  ALONG  THE  2 
COORDINATES,  IN  THE  PRESCRIBED  ORDER.  IF  THE  LITERAL 
CONSTANTS  DESCRIBING  THE  GEOMETRY  ARE  MISS-SPELLED,  AN  ERROR 
MESSAGE  IS  ISSUED,  AND  EXECUTION  STOPPED. 

NOTE  :  THE  2  COORDINATES  ARE  GENERALLY  DENOTED  BY  (X.Y).  IN 
SPHERICAL  R-FYE  GEOMETRY,  FOR  EXAMPLE,  X  MEANS  R,  WHILE 
Y  MEANS  FYE . 

(2)  THE  LEFT,  RIGHT,  BOTTOM,  AND  TOP  BOUNDARIES  ARE  EXTENDED 

1  CELL  BEYOND  THE  LAST  GRID  POINT,  YIELDING  (NX+2) * (NY+2) 
CELLS.  THE  DENSITY  OF  AN  EXTRA  LEFT  CELL  =  LBC  *  (DENSITY 
OF  ADJACENT  CELL  ON  SAME  ROW)  +  RHOLBC.  BY  ADJUSTING  THE 
VALUES  OF  THE  TWO  1-D  REAL  ARRAYS  (  OF  DIMENSION  NY  +  2  )  LBC 
AND  RHOLBC,  VARIOUS  TYPES  OF  BOUNDARIES  CAN  BE  SIMULATED. 
SIMILAR  RELATIONS  APPLY  FOR  RIGHT,  BOTTOM,  AND  TOP 
BOUNDARIES,  DENOTED  BY  R,  B,  AND  T,  RESPECTIVELY.  NOTE 
THAT  BOTTOM  AND  TOP  ARRAYS  ARE  NX+2  CELLS  LONG. 

(3)  ALL  THE  BOUNDARIES  ARE  CONSIDERED  PERMEABLE  TO  DIFFUSION  AND 
ANTI-DIFFUSION  FLUXES,  UNLESS  A  CALL  TO  ENTRY  SOLDFY  INFORMS 
THE  ROUTINE  OTHERWISE.  ANY  OF 

CALL  SOLDFY <  4HLEFT ,  KSTRT,  REND  ' 

CALL  SOLDFY <  4HRITE,  KSTRT,  KEND  ) 

CALL  SOLDFY <  4HB0TM,  KSTRT,  KEND  ) 

CALL  SOLDFY <  3HT0P,  KSTRT,  KEND  ) 

MAKES  THE  LEFT,  RIGHT,  BOTTOM,  OR  TOP  BOUNDARIES  IMPERMEABLE 
TO  BOTH  DIFFUSION  AND  ANTI-DIFFUSION  FLUXES  FROM  CELL 
NUMBER  KSTRT  TO  CELL  NUMBER  KEND,  INCLUSIVE. 

NOTE  :  CELL  1  IS  NOW  THE  EXTRA  CELL  BEYOND  THE  BOUNDARY, 
CONFINING  CELLS  2  TO  NX+1,  OR  NY+I. 

ANY  NUMBER  OF  CALLS  TO  SOLDFY  IS  ALLOWED,  MAKING 
IT  POSSIBLE  TO  SOLIDIFY  UNCONNECTED  PATCHES  ALONG  EACH 
BOUNDARY.  EACH  TIME  SOLDFY  IS  CALLED,  A  MESSAGE  EXPLAINING 
THE  ACTION  TAKEN  IS  ISSUED. 

(4)  CALLS  TO  ENTRY  PRODIC,  FOR  EXAMPLE, 

CALL  PRODIC  (  1  ,  1HX  ) 

INFORM  THE  ROUTINE  TO  TREAT  THE  1  ST  OR  2  ND  COORDINATE  AS 
PERIODIC.  THE  SECOND  ARGUMENT  IS  JUST  TO  GENERATE  A  LABEL > 
THE  MESSAGE  "  X  COORDINATE  FERIODIC  "  IS  ISSUED.  SIMILARLY, 
CALL  PRODIC (  2  ,  3HFYE  ) 

MAKES  THE  2  ND  COORDINATE  PERIODIC,  AND  THE  MESSAGE  "  FYE 
COORDINATE  PERIODIC  "  ISSUED. 

IF  THE  PERIODIC  CALL  IS  MADE  FOR  A  COORDINATE  THAT  SHOULDN'T 
BE  FERIODIC,  A  WARNING  MESSAGE  IS  ISSUED,  THEN  EXECUTION 
PROCEEDS. 
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THE  GRID  IS  INITIALIZED  BY  A  CALL  TO  ENTRY  ORIGRD: 

CALL  ORIGRD <  XGN,  YGN  ) 

WHERE  XGN,  YGN  ARE  TWO  1-D  REAL  ARRAYS  OF  DIMENSIONS  NX+1, 
NY+1,  CONTAINING  THE  LOCATIONS  QF  X ,  Y  INTERFACES. 

ORIGRD  WILL  THEN  CONSIDER  THESE  AS  THE  INITIAL  LOCATIONS. 

AT  THE  BEGINNING  OF  EACH  TIME  STEP 
CALL  NGRID (  XGN,  YGN  ) 

WILL  EVALUATE  VOLUME,  MEAN  INTERFACE  AREA,...  OF  CELLS, 
WHEREAS 

CALL  QGRID  <  XGN,  YGN  ) 

AT  THE  END  OF  EACH  TIME  STEP,  RESET  THE  OLD  ARRAYS 
FOR  THE  NEXT  TIME  STEP. 

IT  IS  ASSUMED  THAT  THE  GRID  IS  MOVING,  UNLESS  CALLS  TO  ENTRY 
FIX GRD  FIX  ONE  OR  BOTH  OF  THE  COORDINATES  GRIDS. 

CALL  FI XGRD <  1  ,  1HX  ) 

INFORMS  THE  ROUTINE  THAT  THE  1  ST  COORDINATE  GRID  IS  FIXED. 
THE  SECOND  ARGUMENT  IS  JUST  TO  GENERATE  A  LABEL;  THE  MESSAGE 
"X  GRID  FIXED  "  IS  ISSUED.  SIMILARLY, 

CALL  FI XGRD (  2  ,  1HZ  ) 

FIXES  THE  2  ND  COORDINATE  GRID  AND  ISSUES  THE  MESSAGE  "  Z 
GRID  FIXED  ".  IF  BOTH  COORDINATES  GRIDS  ARE  FIXED,  CALL 
NGRID,  THEN  QGRID,  ONLY  ONCE  AFTER  INITIALIZATION. 

A  PARTICULAR  ANTI-DIFFUSION  FLUX  CORRECTOR  IS  SELECTED  BY  A 
CALL  TO  ENTRY  SETLMT  : 

CALL  SETLMT <  5HB0RIS,  4HB00K  ) 

INVOKES  BORIS-BOOK  FLUX  LIMITER,  WHILE 
CALL  SETLMT (  7HZALESAK ,  1H  ) 

INVOKES  ZALESAK  FLUX  LIMITER.  THE  ARGUMENTS  REFER  TO  THE 
ORIGINATORS  OF  THE  FLUX  LIMITER.  IF  THE  LITERAL  CONSTANTS 
DESCRIBING  A  LIMITER  ARE  MISS-SPELLED,  AN  ERROR  MESSAGE  IS 
ISSUED,  AND  EXECUTION  STOPPED. 

A  TIME  STEP  STARTS  BY  A  CALL  TO  ENTRY  NGRID,  FOLLOWED  BY 
CALL  VQLFLX  <  U,  V,  DT  > 

WHERE  U,V  ARE  TWO  2-D  REAL  ARRAYS  OF  DIMENSIONS  <NX+2> * (NY+2) 
CONTAINING  THE  COMPONENTS  OF  VELOCITY  VECTOR  AT  THE  CELLS 
CENTERS.  DT  IS  THE  TIME  STEP. 

BEFORE  EACH  CALL  TO  FCT2D,  THE  SOURCE  TERM  IS  DETERMINED 
BY  A  SEQUENCE  OF  CALLS  : 

CALL  CLRSRC 

CLEARS  THE  SOURCE  TERM  WHICH  REMAINS  ZERO  UNTIL  ANY  OF  THE 
NEXT  CALLS  IS  DONE.  EACH  CALL  ADDS  TO  THE  SOURCE  TERM. 

ANY  NUMBER  OF  CALLS  IS  ALLOWED,  TO  FORM  THE  TOTAL  VALUE  OF 
THE  SOURCE  TERM. 

CALL  SOR'CES  (  3HBDF ,  SORCE,  DT  ) 

ADDS  A  BODY  TYPE  FORCE,  WHERE  SORCE  IS  A  2-D  REAL  ARRAY  OF 
DIMENSION  (NX +2) * (NY+2)  CONTAINING  THE  BODY  FORCES  PER  UNIT 
VOLUME. 


CALL  SORCES (  4HXGRD,  SORCE,  DT  ) 

CALL  SORCES <  4HYGRD  ,  SQRCE ,  DT  ) 

ADDS  THE  X  OR  Y  COMPONENTS  OF  THE  GRADIENT  OF  THE  QUANTITY 
IN  ARRAY  SORCE. 

CALL  SORCES (  3HDIV,  SORCE,  DT  ) 

ADDS  THE  DIVERGENCE  OF  THE  QUANTITY  IN  SORCE.  ENTRY  SORCES 
DETERMINES  WHICH  FORM  OF  GRADIENT  OR  DIVERGENCE  TO  USE 
ACCORDING  TO  THE  GEOMETRY.  ALTERNATIVELY,  ONE  CAN 
SEPARATELY  CALL  ENTRY  BODY  FOR  BODY  FORCES,  XGRAD  OR  YGRAD 
FOR  THE  GRADIENT  IN  CARTESIAN  COORDINATES,  RCGRAD  OR  YGRAD 
FOR  THE  GRADIENT  IN  CYLINDRICAL  R-Z  GEOMETRY,...  OR  XGRAD  AND 
YGRAD  FOR  DIVERGENCE  IN  CARTESIAN  GEOMETRY,  RCDIV  AND  YGRAD 
DIVERGENCE  IN  CYLINDRICAL  R-Z  GEOMETRY,... 

(10)  THE  TIME  STEP  ENDS  BY  A  CALL  TO  OGRID 

(11)  FOR  2  ND  ORDER  ACCURACY,  STEPS  (8)  ,  (9)  ARE  PERFORMED  TWICE. 

ONCE  WITH  DT=  TIME  STEP  /  2  FOR  THE  HALF  TIME  STEP,  THEN 

DT  =  TIME  STEP  FDR  THE  WHOLE  TIME  STEP. 

ENTRIES  ... 

ENTRY  NGRID ( XGN, YGN) 

ENTRY  OGRID (XGN, YGN) 

ENTRY  ORIGRD(XGN, YGN) 

ENTRY  VOLFLX (U, V, DT) 

ENTRY  SORCES (SRCTYP, SORCE, DT) 

ENTRY  CLRSRC 
ENTRY  BODY (SORCE, DT) 

ENTRY  XGRAD ( SORCE , DT ) 

ENTRY  YGRAD (SORCE, DT) 

ENTRY  RCGRAD ( SORCE , DT ) 

ENTRY  RCDIV (SORCE, DT) 

ENTRY  SETGOM ( GOMTRY , CRD 1 , CRD2 , N 1 , N2 ) 

ENTRY  PROD I C (CRDNT , CRD) 

ENTRY  SETLMT <  LMTR 1 , LMTR2 ) 

ENTRY  FI XGRD (CRNT , CRD) 

ENTRY  SOLDFY (BONDRY, KSTRT, KEND) 


CALLS  TO  . . . 

SUBROUTINE  NL'MU  (NI  ,  NU,  EPS,  NUV,  MUV) 
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DECLARATIONS  : 


PARAMETER  NPX =100, NPY* 100 

PARAMETER  NP1X=NPX+1 ,NP1Y=NPY+1 

PARAMETER  NP2X=NPX+2 . NP2Y=NPY+2 


INTEGER  GEOM,  CRDNT 
INTEGER  IFLX <NP2X,NP2Y) 


LOGICAL  LSRC 
LOGICAL  SNKRNZ 


LOGICAL 

LOGICAL 


XCHNG i YCHNG 
XPRDC> YPRDC 


C 

c 


REAL  TEXT (20) 

REAL  LMTR1 (2) ,LMTR2(2) 

REAL  TGM (3) ,TCRD<6> ,TLM1 (2,4) ,TLM2(2,4> 

REAL  LBC (NP2Y) , RHOLBC (NP2Y ) 

REAL  RBC (NP2Y) , RHQRBC (NP2Y) 

REAL  BBC ( NP2X ) , RHOBBC ( NP2X ) 

REAL  TBC (NP2X ) , RHQTBC (NP2X ) 

REAL  PRMBLL (NP2Y ) , FRMBLR (NP2Y) , PRMBLB ( NP2X ) , PRMBLT (NP2X ) 

REAL  XGQ(NPIX) , XGN(NRlX) , YGO(NPIY) , YGN(NPIY) 

REAL  XG(NPIX) , DXG(NPIX) , YC(NPIY) , DYG(NPIY) 

REAL  DXGO (NP2X ) , DXGN(NP2X) , DYGO (NP2Y) , DYGN (NP2Y) 

REAL  RDXGN(NP2X) ,RDYGN(NP2Y) 

REAL  DXGNH(NPIX)  ,  RDXGNH  (NP1 X )  ,  DYGNH  (NF'IY)  ,  RDYGNH  (NP1Y) 

REAL  AX (NPX) , AY(NPY) 

REAL  SQ(NPIX) ,S00(NP1X> ,S0N(NP1X) 

REAL  RHO ( NP2X , NP2Y , 2 ) 

REAL  U (NP2X , NP2Y ) , ADUDT (NP1 X , NPY ) 

REAL  V (NP2X , NP2Y) , ADVDT (NPX , NPlY) 

REAL  CELMAS (NPX , NPY) , SOURCE (NP2X , NP2Y) , SORCE (NP2X , NP2Y ) 

REAL  TEMPI (NP2X, NP2Y) , TEMP2 (NP2X , NP2Y) 

REAL  TEMP3 (NP1 X , NP1Y) , TEMP4 (NP1X , NPlY) 

REAL  TEMPS  (NP2X,NF'2Y)  ,  TEMP6  (NP2X  ,  NF'2Y) 

REAL  OLD  VOL  ( NF'2X  »  NP2Y )  ,  RVOL  (NP2X  ,  NP2Y) 

REAL  A  V  X  VL  ( NP 1  X  ,  NP  1 Y ),  F:  A  V  X  VL  <  NP  1  X  ,  NP  1 Y  > 

REAL  AVYVL  ( NP  1  X ,  NP  1 Y )  ,  F.AVYVL  ( NP  1  X  ,  NP  1  Y  > 

REAL  XMSFLX  (NP1X.NP1Y)  ,  YMSFLX  (NP  IX,  NF'IY) 

REAL  XDFFLX (NP2X,NP2Y) , YDFFLX (NP2X,NP2Y) 

REAL  XNTFLX  (NP1  X  ,  NP1 Y)  ,  YNTFLX  (NF'IX,  NPlY) 
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REAL  EF'SX  (NP1X.NP1Y)  ,  NUX  (NP2X,NP2Y)  ,  MUX  (NP2X.NP2Y) 

REAL  EPS Y  ( NP 1 X , NP 1 Y )  , NUY (NP2X , NP2Y) , MUY (NP2X , NP2Y ) 

REAL  NUX  VOL  ( NP2X  >  NF'2Y )  ,  MUX  VOL  (NP2X ,  NP2Y) 

REAL  NUYVOL ( NP2X  ,  NP2Y ) ,  MUYVOL  (NP2X  ,  NF'2Y> 

C 

REAL  MXFLX (NP2X, NF'2Y) , MNFLX (NP2X, NP2Y) 

REAL  FLX IN (NP2X  ,  NP2Y) . FLXOUT ( NP2X  »  NP2Y ) 

REAL  RHOMX (NP2X,NP2Y) , RHOMN (NP2X , NP2Y) 

REAL  MX  IN (NP2X . NP2Y)  , MXOUT (NP2X , NP2Y) 

REAL  DIFF  (NP2X  ,  NP2Y )  ,  FLX  ( NF'2X  ,  NP2Y ) 

C 

REAL  R I N  <  NF'2X  ,  NP2Y )  ,  ROUT  ( NP2X  ,  NF'2 Y ) 

REAL  XFLXCR (NP2X, NP2Y) , YFLXCR (NP2X  ,  NP2Y ) 


EQUIVALENCE 

EQUIVALENCE 

EQUIVALENCE 

EQUIVALENCE 

EQUIVALENCE 

EQUIVALENCE 

EQUIVALENCE 

EQUIVALENCE 

EQUIVALENCE 


(TEMPI, FLXIN,RIN> 

( TEMP2 , FLXOUT , ROUT ) 

(TEMPS,  XMSFLX ,  XNTFLX  ,  AVXVL ,  RAVXVL,  EF'SX ) 

( TEMP4 ,  YMSFLX ,  YNTFLX  ,  AVYVL ,  RAVYVL ,  EF’SY ) 
(TEMPS, XDFFLX , YDFFLX , NUX , NUXVOL , NUY , NUYVOL ) 
( TEMP5, OLDVOL , RVOL , DIFF , MXFLX , FLX , IFLX ) 

( TEMP6 , SOURCE , MUX , MUX VOL , MUY , MUYVOL ) 

( TEMP6 , RHOMX , RHQMN , M  X I N , M  X  OUT , MNFL X ) 

(TEMP 6, XFLXCR, YFLXCR) 


DATA  TEXT ( 1 ) , TEXT (2) , TEXT (3) /4H  M, 4HISS-, 4HSPEL/ 
DATA  TEXT (4) , TEXT (5) /4HLING , 4H  OF  / 

DATA  TEXT (9) »  TEXT ( 10)  , TEXT (11) /4H  IDE, 4HNTIF, 4HIER  / 
DATA  TEXT ( 12) , TEXT ( 13) , TEXT ( 14) /4H  GE , 4H0MET , 4HRY  / 
DATA  TEXT(15) , TEXT (1 6 ) , TEXT (1 7) /4HFLUX , 4H  L IM, 4HI TER/ 
DATA  TEXT (18) , TEXT (19) , TEXT (20) /4HS0UR, 4HCE  T,4HYFE  / 

DATA  TOM ( 1 ) , TGM (2) , TGM ( 3) /4HCART , 4HCYL  , 4HSPH  / 

DATA  TCRD ( 1 ) , TCRD (2) , TCRD (3) /4HX  , 4HY  , 4HZ  / 
DATA  TCRD (4) , TCRD (5) , TCRD (6) /4HR  , 4HCETA , 4HFYE  / 

DATA  TBND ( 1 ) , TBND (2) /4HLEFT , 4HRITE/ 

DATA  TBND (3) , TBND (4) /4HB0TM, 4HT0P  / 

DATA  LMTR1 , LMTR2/4*4H  / 

DATA  TLM1 (1,1), TLM1 (2,1) /4HB0RI , 4HS  / 

DATA  TLM2 (1,1), TLM2 (2,1) /4HB00K , 4H  / 

DATA  TLM1 (1,2) , TLM1 (2,2) /4HZALE, 4HSAK  / 

DATA  TLM2 (1,2) , TLM2 (2, 2) /4H  , 4H  / 

DATA  BDF/3HBDF/, XGRD, YGRD/4HXGRD, 4HYGRD/ >  DIV/3HDIV/ 
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FORMATS  : 


10  FORMAT (///5X,7HWARNING,5X, 11A4) 

20  FORMAT  ( ///5X  >  7HWARNING,  5X »  A4 ,  2X  *  22HSH0ULD  NOT  BE  PERIODIC) 

30  FORMAT ( ///5X , 24HALL  BOUNDARIES  PERMEABLE) 

40  FORMAT ( ///5X, A4, 2X . 19HC00RD I NATE  PERIODIC) 

50  FORMAT (///5X.A4.2X, 10HGRID  FIXED) 

60  FORMAT ( / / /5X  »  A4 , 2X , 22HB0UNDARY  SOLID  BETWEEN, 

&  2X , 4HCELL >  2X ,  14, 2X , 3HAND, 2X , 14) 

70  FORMAT < / / /5X , 25HGE0METRY  NOT  INCLUDED  YET) 
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EVALUATE  OLD  CELL  MASS  "CELMAS" 

DO  110  J= 1 »  NY 
DO  110  1=1, NX 

OLDVOL ( I , J) =DXG0 ( 1+1 ) *DYG0  < J+l ) 

CELMAS  <  I ,  J )  =RHQ  ( 1  +  1 ,  J+l ,  1  >  *0LDV0L  ( I ,  J ) 

ADD  SOURCE  TERM  "SOURCE'*  WHEN  APPROPRIATE 
IF  ( . NOT. LSRC)  GO  TO  125 


DO  120  J=1,NY 
DO  120  1=1, NX 

CELMAS ( I , J) =CELMAS < I , J) -SOURCE ( I ,  J) 
CONTINUE 


EVALUATE  X-CONVECTION  FLUX  "XMSFLX" 

DO  130  J=1 , NY 
DO  130  I =1 , NXP1 

TEMPS ( I , J) =RHO ( 1+1 , J+l , KO) +RHO ( I , J+l , KO) 
TEMP5 ( I , J) =0. 5#ADUDT ( I , J) 

XMSFLX ( I , J) =TEMP3 ( I , J) * TEMPS ( I , J) 


EVALUATE  Y-CONVECTION  FLUX  " YMSFLX " 

DQ  140  J= 1 , NYP 1 
DO  140  1=1, NX 

TEMP4 < I , J) =R HO ( I  + 1 , J+ 1 , KO) +RHO ( I  + 1 , J , KO  > 

TEMP6 ( I , J) =0. 5*ADVDT ( I , J) 

140  YMSFLX (I, J)=TEMP4(I, J) *TEMP6(I, J) 

EVALUATE  X  AND  Y  TRANSPORTED  DENSITIES 
DO  150  J= 1 , NY 
DO  150  1=1, NX 

TEMP5 ( I , J) =XMSFLX (I, J) -XMSFLX ( 1+1 , J) 

TEMP6 < I , J) =YMSFLX ( I , J) -YMSFLX ( I , J+l ) 

TEMPS ( I , J) =C£LMAS < I , J) +TEMP5 ( I , J) 

TEMP4 ( I , J) =CELMAS ( I , J) +TEMP6 (I , J) 

CELMAS  < I , J) =TEMP3 ( I , J) +TEMP6 ( I , J) 

RVOL ( I , J) =RDXGN ( 1+1 >  *RDYGN ( J+l ) 

RHO ( 1+1 , J+l , KN) =TEMP3 ( I , J) *RVOL ( I , J) 

150  TEMP4  v  I ,  J )  =TEMF'4  ( I ,  J )  *RVOL(I, J) 

EVALUATE  X-TRANSPORTED  DENSITY  AT  LEFT  AND  RIGHT  BOUNDARIES 
DO  170  J=2 , NYP 1 

RHO ( 1 , J, KN) =LBC ( J ) *RHQ ( IL , J , KN ) +RHOLBC < J ) 

170  RHO (NXP2, J, KN) =RBC< J) *RHO  < IR, J, KN) +RHORBC ( J ) 
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C  EVALUATE  "EF'SX" 

DO  130  J= 1 i NY 
DO  180  1  =  1,  NXF‘1 

RAVXVL ( I , J) =RDXGNH ( I ) #RDYGN ( J+l > 

180  EPSX ( I ,  J) =RAVXVL (1 , J) *ADUDT <  I ,  J) 

EVALUATE  X  DIFFUSION  AND  ANT  I -D I FFUS I ON  TEOEFF I C I ENTS  "NUX",  "MUX" 
CALL  NUMU (NXP1 , NY , EPSX  >  NUX  »  MUX ) 

CANCEL  THE  DIFFUSION  AND  ANTI-DIFFUSION  X-FLUXES  THROUGH  SOLID 
PORTIONS  OF  LEFT  AND  RIGHT  BOUNDARIES 
DO  185  J=1,NY 

NUX ( 1 , J) =NUX <  1 ,  J) *PRMBLL <  J+l  > 

MUX ( 1 >  J) =MUX ( 1 1 J) *PRMBLL(J+1) 

NUX (NXP1 , J) =NUX (NXP1, J) *PRMBLR ( J+ 1 ) 

185  MUX (NXP1 *  J) =MUX (NXPl , J) *PRMBLR (J+l ) 

EVALUATE  X  DIFFUSION  AND  ANTI-DIFFUSION  FLUXES  "XDFFLX"  ,  ” XNTFLX " 

DO  190  J=1 >  NY 
DO  190  1=1, NXPl 
AVXVL  ( I ,  J)  =DXGNH  <  I )  *DYGN  <  J+l ) 

NUXVOL (1 i J) =NUX (I, J) f AVXVL ( I , J ) 

MUXVOL ( I , J) =MUX ( I , J) *AVXVL ( I , J) 

TEMPS ( I , J) =RHO ( 1+1 , J+l , KO) -RHO ( I , J+l , KO) 

XDFFLX  ( I ,  J)  =NUXVOL  ( I ,  J)  *TEMF'3  ( I ,  J) 

TEMPS  <  I ,  J)  =RHQ  ( 1  +  1 ,  J+l ,  KN)  -RHO  ( I  ,  J+l ,  K.N> 

190  XNTFLX ( I , J) =MUXVOL  < I , J) * TEMPS ( I , J) 

ADD  X-DIFFUSION  TO  "CELMAS" 

DO  200  J= 1 , NY 
DO  200  1=1, NX 

RHO <1+1 , J+l ,KN)=TEMP4(I, J) 

TEMP6 ( I , J) = XDFFLX ( 1+1 , J) -XDFFLX ( I , J) 

200  CELMAS  < I , J) =CE LMAS (I,J) +TEMP6 (I , J> 

EVALUATE  Y-TRANSPORTED  DENSITY  AT  BOTTOM  AND  TOP  BOUNDARIES 
DO  210  1=1 , NXP2 

RHO (1,1, KN) =BBC ( I ) *RHO ( I , JB, KN) +RHOBBC ( I ) 

210  RHO ( I , NYP2, KN) =TBC ( I ) *RHQ ( I , JT, KN) +RHOTBC ( I ) 

EVALUATE  "EF'SY" 

DO  220  J=  1 ,  NYF'  1 
DO  220  1=1, NX 

RAVYVL ( I , J) =RD XGN < 1+ 1 ) *RDYCNH(J) 

220  EF’SY  ( I ,  J)  =P.AVYVL  ( I ,  J)  *ADVDT  ( I ,  J) 

EVALUATE  Y  DIFFUSION  AND  ANTI-DIFFUSION  COEFFICIENTS  "NUY" ,  "MUY" 
CALL  NUMU  < NX , NYP 1 , EPSY , NUY , MUY ) 
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CANCEL  THL  DIFFUSION  AND  ANTI-DIFFUSION  Y-FLUXES  THROUGH  SOLID 
PORTIONS  OF  BOTTOM  AND  TOP  BOUNDARIES 
DO  225  1=1, NX 

NUY (1,1) =NUY (1,1) *PRMBLB <  I  + 1  ) 

MUY (1,1) =MUY (1,1) *PRMBLB <  I  + 1  > 

NUY ( I , NYP1 ) =NUY ( I , NYP1 )  *P RMBLT ( 1  +  1 ) 

225  MUY ( I , NYP1 ) =MUY ( I , NYP1 ) *PRMBLT (1  +  1  ) 

EVALUATE  Y  DIFFUSION  AND  ANTI -DIFFUSION  FLUXES  " YDFFLX "  ,  " YNTFLX " 

DO  230  J=1 , NYP1 
DO  230  1=1, NX 

AVYVL ( I , J) =DXGN  < 1  +  1 ) *DYGNH  < J) 

NUYVOL ( I , J) =NUY ( I , J) *AVYVL  < I , J) 

MUYVOL ( I , J) =MUY ( I , J) *AVYVL ( I , J) 

TEMP4 ( I , J) =RH0 ( I  + 1 , J+ 1 , KO  > -RHO ( I + 1 , J , KO ) 

YDFFLX ( I , J) =NUYV0L ( I , J) *TEMP4 ( I , J) 

TEMP4 ( I , J) =RHO ( 1  +  1 >  J+l , KN) -RHO ( I + 1 , J , KN ) 

230  YNTFLX ( I , J) =MUYVQL ( 1 , J) *TEMP4 ( I ,  J) 

ADD  Y-DIFFUSION  TO  “CELMAS" 

DO  240  J= 1 , NY 
DO  240  1=1, NX 

TEMP6U,  J)=YDFFLX  ( I ,  J+ 1 ) -YDFFLX  (I,  J) 

240  CELMAS ( I ,  J ) =CELMAS < I , J ) +TEMP6 ( I , J  > 

IF  SYNCHRONIZATION  OF  ANTI-DIFFUSION  FLUXES  IS  SPECIFIED,  SKIP 
EVALUATION  OF  CORRECTION  FACTORS 
IF ( SNKRNZ )  GO  TO  445 


EVALUATE  TRANSPORTED-DIFFUSED  DENSITY 
DO  250  J= 1 , NY 
DO  250  1=1, NX 

RVOL ( I , J ) =RDXGN < 1+ 1 ) CRDYGN <  J+ 1  > 

250  RHO (1+1 , J+l , KN) =CELMAS < I , J) *RVGL ( I , J) 

EVALUATE  TRANSPORTED-DIFFUSED  DENSITY  AT  BOTTOM  AND  TOP  BOUNDARIES 
DO  260  1=2, NXP1 
TEMPI (I,  1)=1.0 
TEMP2 ( I - 1 , NYP 1 ) = 1 . 0 

RHO  <  1 ,  1 ,  KN)  =BBC  ( I )  *RHO  (I ,  JB,  KN)  +RHDBBC  ( I ) 

260  RHO ( I » NYP2, KN) =TBC ( I ) *RHO  v I , JT , KN) +RHOTBC ( I ) 

EVALUATE  TRANSPORTED-DIFFUSED  DENSITY  AT  LEFT  AND  RIGHT  BOUNDARIES 
DO  270  J=2  >  NYP1 
TEMPI (1, J)=1.0 
TEMP2  (NXF'l  ,  J-l )  =1 . 0 

RHO (1 , J,KN) =LBC ( J) IRHO ( IL, J, KN) +RHOLBC ( J) 

270  RHO (NXP2, J, KN) =RBC ( J) *RHQ ( IR, J,KN) +RHORBC ( J) 


A10 


nnnnn  nnnn  n  n  nn 


CANCEL  THE  ANTI-DIFFUSION  X-FLUX  IF  IT  IS  OPPOSITE  TO  ITS  LOCAL 
TRANSPORTED-DIFFUSED  DENSITY  GRADIENT  AND  ANY  OF  THE  ADJACENT  ONES 
DO  280  J=1,NY 
DO  280  1  =  1 >  NXP 1 

DIFF  (  I , J) =RHO ( I  + 1  >  J+ 1 1 KN) -RHO ( I , J+l ,  KN) 

DO  290  J=1,NY 
DO  290  1=1, NX 

TEMPI  ( 1  +  1 ,  J+l )  =XOR  ( XN  IFLX  ( I-f-1 ,  J)  ,  DIFF  ( I ,  J)  ) 

TEMP2CI, J) =XGR (XNTFLX (I, J) , DIFF < 1+1 , J) ) 

DO  300  J=1 , NY 
DO  300  1=1 ,NXP1 

TEMP5<I, J)=XDR< XNTFLX (I, J> , DIFF (I, J) ) 

TEMP6 ( I , J) =QR (TEMPI (I, J+l) , TEMP2 ( I , J) ) 

FLX ( I , J) =AND (TEMP5  < I , J) , TEMP6 ( I , J ) ) 

FLX ( I , J) =COMPL (FLX ( I , J) ) 

I FLX (I, J)=LSHF(IFLX(I, J) ,-31) 

FLX ( I , J) =FLQAT ( IFLX ( I , J) ) 

XNTFLX (I, J) =XNTFLX (I, J) *FLX  (I,  J) 

IF  X— COORD I NATE  IS  PERIODIC  AND  EITHER  LEFT  OR  RIGHT  BOUNDARY'S 
ANTI-DIFFUSION  FLUX  IS  CANCELLED,  CANCEL  THE  OTHER 
IF (. NOT. XPRDC)  GO  TO  305 

DO  304  J= 1 , NY 

XNTFLX (1, J)=AND( XNTFLX (1, J) , XNTFLX (NXP1 , J) ) 

XNTFLX (NXP 1 , J) =XNTFLX ( 1 , J) 

CONTINUE 


CANCEL  THE  ANT I -DIFFUSION  Y-FLUX  IF  IT  IS  C 
TRANSPORTED-DIFFUSED  DENSITY  GRADIENT  AND  P 
DO  310  J= 1 , NYP 1 
DO  310  1=1, NX 

D I FF  < I , J ) =RHO ( I  + 1 , J+ 1 , KN ) -RHO ( I  + 1 , J , KN ) 


DO  320  J= 1 , NY 
DO  320  1=1, NX 

TEMPI (1+1, J+1)=XQR(YNTFLX(I, J+l) , DIFF (I, J) ) 
TEMP2 ( I , J) =XOR ( YNTFLX (I, J) , DIFF (I, J+l) ) 

DO  330  J=1 , NYP 1 
DO  330  1=1, NX 

TEMPS  < I , J) =XOR (YNTFLX ( I , J) , DIFF  < I , J) ) 

TEMP6 ( I , J) =0R ( TEMP 1 < I + 1 , J ) , TEMP2 ( I , J) ) 

FLX ( I , J) =AND (TEMPS ( I , J) , TEMP6 ( I , J ) ) 

FLX (I, J)=COMPL(FLX (I, J) ) 

IFLX ( I , J) =LSHF ( IFLX ( I , J) ,-31) 

FLX ( I , J) =FLOAT ( IFLX ( I , J) ) 

YNTFLX  < I , J) =YNTFLX ( I , J) *FLX ( I , J) 


OPPOSITE  TO  ITS  LOCAL 
ANY  OF  THE  ADJACENT  ONES 


C  IF  Y-COORDINATE  IS  PERIODIC  AND  EITHER  BOTTOM  OR  TOP  BOUNDARY’S 

C  ANTI-DIFFUSION  FLUX  IS  CANCELLED,  CANCEL  THE  OTHER 

IF  ( .  NOT.  YF'RDC)  GO  TO  335 
C 

DO  334  1=1, NX 

YNTFLX (1,1) =AND ( YNTFLX (1,1), YNTFLX ( I , NYP1 ) ) 

334  YNTFLX ( I , NYP1 ) =YNTFLX < I ,  1  ) 

C 

335  CONTINUE 
C 

C 

C  EVALUATE  NET  INCOMING  "FLX IN" ,  OUTGOING  "FLXQUT"  ANTI-DIFFUSION 

DO  340  J= 1 , NY 
DO  340  1=1, NXP 1 

TEMP5 ( I , J) =ASHF ( XNTFLX ( I , J) , -31 ) 

MNFLX (I , J) =AND( XNTFLX ( I, J) , TEMPS ( I , J) ) 

TEMPS ( I, J)=XOR( XNTFLX (I, J) , TEMP5 ( I , J) > 

340  MXFLX ( I , J)=AND (XNTFLX ( I ,  J) , TEMPS ( I ,  J) ) 

C 

DO  350  J= 1 , NY 
DO  350  1=1, NX 

FLX IN ( 1+1 , J+l ) =1 . E-50+MXFLX ( I , J) 

FLXOUT ( 1  +  1 , J+l ) =1 . E-50-MNFLX ( I , J) 

FLX IN ( 1  +  1 , J+l ) =FLX IN  < 1  +  1 , J+l) -MNFLX ( 1  +  1 , J) 

350  FLXQUT ( 1+1 , J+l ) =FLXOUT ( 1+1 , J+l ) +MXFLX ( 1+1 , J) 

C 

DO  360  J= 1 , NYP 1 
DO  360  1=1, NX 

TEMPS ( I , J) =ASHF (YNTFLX ( I , J)  ,  -31 ) 

MNFLX ( I , J) =AND ( YNTFLX ( I , J )  ,  TEMPS ( I , J ) ) 

TEMP5(I, J)=XOR< YNTFLX (I, J) , TEMPS ( I ,  J) ) 

360  MXFLX (I, J)=AND (YNTFLX (I, J) , TEMP5 ( I , J ) ) 

C 

DO  370  J= 1 , NY 
DO  370  1=1, NX 

FLX IN ( I +  1 , J+l ) =FLX IN ( 1  +  1 , J+l >  +MXFLX ( I, J) 

FLXOUT (1+1 , J+l )=FLXOUT ( 1+1 , J+l ) -MNFLX ( I , J) 

FLX I N ( I + 1 , J+ 1 ) =FLX IN ( 1+1 , J+l ) -MNFLX ( I , J+ 1 ) 

370  FLXOUT ( 1+1 , J+l ) =FLXOUT ( 1+1 , J+l ) +MXFLX ( I , J+l ) 

C 

C 

GO  TO  (375,385)  ILMTR 
C 

37S  CONTINUE 

C  IF  BORIS-BOOK  FLUX  LIMITER  IS  REQUESTED,  USE  TRANSPORTED-DIFFUSED 

C  DENSITY  TO  BOUND  NEW  DENSITY 

DO  380  J= 1 , NYP2 
DO  380  1  =  1 ,  NXF'2 
380  TEMPS ( I , J) =RHQ ( I , J, KN) 

C 

GO  TO  395 
C 


A12 


S3 


385  CONTINUE 

C  IF  ZALESAK  FLUX  LIMITER  IS  REQUESTED,  USE  MAXIMUM  OF  OLD  AND 

C  TRANSPORTED-DIFFUSED  DENSITY  AS  UPPER  BOUND  FOR  NEW  DENSITY 

DO  390  J=1 , NYP2 
DO  390  1  =  1 , NXP2 

390  TEMPS  ( I ,  J)  =AMAX  1  (RHO  ( I ,  J ,  KO )  ,  RHO  ( I ,  J.KN)  ) 

C 

395  CONTINUE 

C 

DO  400  J=2>  NYF‘1 

C  EVALUATE  MAXIMUM  ADMISSIBLE  ANTI-DIFFUSION  INTO  CELL  "MXIN",  AND 

C  IN  TURN  CORRECTION  FACTOR  "RIN" 

DO  400  1=2, NXP1 

TEMP6 ( I , J) =AMAX 1 (TEMPS ( I, U) , TEMP5 ( 1-1 , J) ) 

TEMP6 ( I , J) =AMAX 1 (TEMF'6 ( I ,  J)  ,  TEMPS ( 1  +  1 ,  J)  ) 

TEMF'6 ( I , J ) =AMAX 1 (TEMF'6  ( I,  J) , TEMPS ( I , J- 1 ) ) 

RHOMX ( I , J) =AMAX 1 (TEMP6(I, J) , TEMPS ( I , J+l ) ) 

TEMP6 ( I , J) =RHOMX ( I , J) -RHO ( I , J, KN> 

TEMP6 ( I , J) =TEMP6 (I , J) *DXGN ( I ) 

MX  I N ( I , J ) =TEMP6 ( I , J ) *DYGN(J) 

TEMP6 ( I , J) =MX I N ( I , J) /FLXIN ( I , J) 

400  RIN ( I , J) = AM INI (1.0, TEMP6 ( I , J) ) 

C 

GO  TO  (415,405)  ILMTR 
C 

405  CONTINUE 

C  IF  ZALESAK  FLUX  LIMITER  IS  REQUESTED,  USE  MINIMUM  OF  OLD  AND 

C  TRANSPORTED-DIFFUSED  DENSITY  AS  LOWER  BOUND  FOR  NEW  DENSITY 

DO  410  J=1 , NYP2 
DO  410  1  =  1  *  NXP2 

410  TEMF'S  <  I ,  J)  =AMIN1  (RHO  ( I ,  J,  K0>  .RHO  (I,  J.KN)  > 

C 

415  CONTINUE 

C  EVALUATE  MAXIMUM  ADMISSIBLE  ANTI-DIFFUSION  OUT  OF  CELL  "MXOUT", 

C  AND  IN  TURN  CORRECTION  FACTOR  "ROUT" 

DO  420  J=2 , NYP 1 
DO  420  I =2 , NXP 1 

TEMF'6  ( I ,  J )  =AM  INI  (TEMPS  ( I,  J)  ,  TEMPS  ( I- 1 ,  J )  ) 

TEMP6 ( I , J ) =AM INI ( TEMP6 ( I , J ) , TEMPS ( 1+1 , J) ) 

TEMP6  ( I ,  J )  =AM  INI  (TEMF'6  ( I,  J)  ,  TEMPS  ( I ,  J- 1 )  ) 

RHOMN ( I , J ) =AM I N 1 ( TEMP6 ( I , J) , TEMPS ( I , J+ 1 ) ) 

TEMP6 ( I ,  J) =RHO  < I , J, KN) -RHOMN ( I , J) 

TEMF'6  ( I ,  J)  =TEMP6  ( I ,  J)  *DXGN  ( I ) 

MXOUT ( I , J) =TEMP6 ( I , J) *DYGN ( J) 

TEMP6 ( I , J) =MXQUT ( I , J) /FLXOUT ( I , J) 

420  ROUT  (I,  J)=AMIN1  ( 1 . 0,  TEMF'6  ( I ,  J )  ) 

C 
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IF  A  COORDINATE  IS  NOT  PERIODIC,  "RIN",  ROUT  "  ARE  ASSUMED  TO  BE 
CONTINUOUS  THROUGH  ITS  NORMAL  BOUNDARY 
DO  430  1=2, NXP1 
RIN ( I , 1 ) =RIN ( I , JB) 

RIN ( I , NYP2) =RIN ( I , JT) 

ROUT (I, 1)=R0UT(I, JB) 

430  ROUT ( I , NYP2) =ROUT  < I , JT) 

DO  440  J=2, NYP1 
RIN ( 1 , J) =RIN ( IL, J) 

RIN (NXP2, J) =RIN ( IR, J) 

ROUT ( 1 , J) =ROUT  < IL, J) 

440  ROUT (NXP2, J) =RQUT ( IR, J) 

445  CONTINUE 


LIMIT  ANTI-DIFFUSION  FLUXES  USING  MINIMUM  OF  ADJACENT  CELLS' 
MAXIMUM  ADMISSIBLE  FLUXES 
DO  450  J=1 , NY 
DO  450  1=1 , NXP1 
FLX (I, J)=XNTFLX (I , J) 

IFLX (I, J)=LSHF(IFLX (I, J) ,-31) 

FLX ( I , J) =FLOAT ( IFLX  < I , J) ) 

RHQ ( I  + 1 »  J+ 1 , KN ) =AM I N 1 (RIN (I, J+l) , ROUT < 1+1 , J+l ) ) 

XFLXCR ( I , J) =FLX ( I , J) *RH0(I+1, J+l , KN) 

RHO ( 1+1 , J+l , KN) =AMIN1 (RIN ( 1+1 , J+l ) .ROUT (I, J+l) > 

FLX (I, J)=1.0-FLX(I, J) 

FLX ( I , J) =FLX ( I , J) *RHO ( 1+1 , J+l , KN) 

XFLXCR  < I , J) =XFLXCR ( I , J) +FLX ( I , J) 

450  XNTFLX < I , J) =XNTFLX (I , J) *XFLXCR < I , J) 

DO  460  J= 1 , NYP 1 

DO  460  1=1, NX 

FLX (I, J)=YNTFLX (I, J) 

IFLX (I , J) =LSHF < IFLX (I , J) , -31 ) 

FLX ( I , J) =FLOAT ( IFLX ( I , J) ) 

RHO ( 1+1 , J+l , KN) =AMINl (RIN(I+1,J> , ROUT ( I +1 , J+l ) ) 

YFLXCR ( I , J) =FLX ( I , J) *RHO ( 1+1 , J+l , KN) 

RHO ( 1  +  1 , J+l ,  KN) =AMIN1 (RIM ( 1  +  1 , J+l ) , ROUT ( I  + 1 , J ) ) 

FLX ( I , J) =1 . O-FLX ( I , J) 

FLX  (I ,  J)=FLX  <  I,  J)  '4!  RHO  ( 1  +  1 ,  J+l ,  KN) 

YFLXCR ( I , J) =YFLXCR ( I , J) +FLX ( I , J) 

460  YNTFLX ( I , J) =YNTFLX < I , J) * YFLXCR ( I , J) 
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C  ADD  CORRECTED  ANTI-DIFFUSION  FLUXES  AND  EVALUATE  NEW  DENSITY 
DO  470  J=1 ,  NY 
DO  470  1=1, NX 

TEMP5 ( I , J) =XNTFLX ( I , J) -XNTFLX ( 1  +  1 >  J) 

TEMP6 ( I , J) =YNTFLX ( I , J) -YNTFLX < I , J+l ) 

CELMAS ( I ,  J) =CELMAS ( I , J) +TEMP5 ( I , J) 

CELMAS ( I , J) =CELMAS ( I , J) +TEMP6 ( I , J) 

RVCL (1,J) =RDXGN < 1+1 ) *RDYGN ( J+ 1 ) 

470  RHO( 1+1 , J+l ,KR) =CELMAS (I, J) *RVOL (I , J) 

EVALUATE  NEW  DENSITY  AT  BOTTOM  AND  TOP  BOUNDARIES 
DO  490  1=2, NXP1 

RHO<I, 1 ,  KR) =BBC ( I ) *RHO ( I , JB, KR) +RHOBBC ( I ) 

490  RHO ( I , NYP2, KR) =TBC  < I ) #RHO < I , JT, KR) +RHOTBC ( I ) 

EVALUATE  NEW  DENSITY  AT  LEFT  AND  RIGHT  BOUNDARIES 
DO  500  J=2, NYP1 

RHO  < 1 , J, KR) =LBC ( J) *RHO ( IL , J, KR) +RHOLBC ( J) 

500  RHO (NXP2, J, KR) =RBC ( J) *RHO ( IR, J, KR) +RHORBC ( J) 

RETURN 
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ENTRY  NGRID ( XGN, YGN) 


EVALUATE  AVERAGE  (  BETWEEN  OLD  AND  NEW  )  INTERFACE  VELOCITY  AND  AREf 
AND  NEW  AND  AVERAGE  VOLUME  COMPONENTS. 

IF  X-GRID  OR  Y-GRID  IS  NOT  MOVING.  USE  ITS  OLD  VALUES. 

INTERFACE  VOLUME  IS  CONSIDERED  AVERAGE  OF  ADJACENT  CELLS'  VOLUMES. 


GO  TO  (510.520.530,540,550,560,570)  GEOM 


510  CONTINUE 
CARTESIAN  COORDINATES 

IF ( . NOT. XCHNG)  GO  TO  513 

DO  511  1  =  1  *  NXP1 
XG ( I ) =0. 5* (XGN ( I ) +XGO ( I > ) 

511  DXG(I)=XGN<I>-XGO<I) 

DO  512  1=2, NXP1 
DXGN ( I ) =XGN ( I ) -XGN ( I -1 ) 

AX (I-1)=XG(I)-XG(I-1) 

512  RDXGN ( I ) =1 . O/DXGN ( I ) 

513  CONTINUE 
IF  < . NOT. YCHNG)  GO  TO  580 

DO  516  J=  1 ,  NYF'  1 
YG(J)=0.5*<YGN<J)+YG0(J) ) 

516  DYG ( J) =YGN ( J) — YGO ( J) 

DO  517  J=2,  NYF'l 
DYGN ( J) =YGN ( J ) -YGN  < J- 1 ) 

AY ( J-l ) =YG ( J) -YG ( J-l ) 

517  RDYGN ( J) =1 . O/DYGN ( J) 

GO  TO  580 

520  CONTINUE 
CYLINDRICAL  R-Z  COORDINATES 

IF (.NOT. XCHNG)  00  TO  523 

DO  521  1=1, NX PI 
SCO ( I ) =0.5* (XGO(I) *XGO ( I) ) 
SON ( I ) =0. 5* ( XGN ( I ) *  XGN ( I ) ) 
XG(I)=SaN(I)+SOO(I) 

DXG ( I ) =SQN  < I ) -SQO ( I ) 

521  SQ ( I ) =SQRT ( XG ( I ) ) 

C 
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DO  522  1=2, NXP 1 

DXGN ( I ) =SQN ( I ) -SON  < I - 1 ) 

AX (I-l>=0.5* (XG(I)-XG(I-l) ) 

522  RDXGN ( I ) =1 . Q/DXGN ( I ) 

C 

523  CONTINUE 

IF (.NOT. VCHNG)  GO  TO  580 
C 

DO  526  J=1 , NYP1 
YG<J)=0.5*<YGN(J)+YG0<J>  > 

526  DYG  < J) =YGN  <  J) -YGO  <  J) 

C 

DO  527  J=2, NYP1 

DYGN ( J) =YGN ( J) -YGN ( J-l ) 

AY ( J-l ) =YG ( J) -YG  < J-l ) 

527  RDYGN ( J) =1 . O/DYGN  < J) 

C 

GO  TO  580 
C 

530  CONTINUE 

C  CYLINDRICAL  R-FYE  COORDINATES 

540  CONTINUE 

C  CYLINDRICAL  Z-FYE  COORDINATES 

550  CONTINUE 

C  SPHERICAL  R-THETA  COORDINATES 

560  CONTINUE 

C  SPHERICAL  R-FYE  COORDINATES 

570  CONTINUE 

SPHERICAL  THETA-FYE  COORDINATES 

PRINT  70 
STOP 
C 
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580  CONTINUE 

IF (.NOT. XCHNG)  GO  TO  586 
C 

DXGN ( 1 ) =DXGN ( IL) 

DXGN (NXP2) =DXGN ( IR) 

RDXGN  < 1 ) =RDXGN ( IL) 

RDXGN (NXP2) -RDXGN ( IR) 

C 

DO  585  1  =  1 j  NXP1 

DXGNH  <  I )  =0. 5-K  (DXGN  ( I )  +DXGN  ( 1  +  1 )  ) 

585  RDXGNH ( I ) =0. 5# (RDXGN ( I ) +RDXGN ( 1+1 ) ) 
C 

586  CONTINUE 

IF  .'.NOT.  YCHNG)  RETURN 
C 

DYGN ( 1 ) =DYGN ( JB ) 

DYGN (NYP2) =DYGN ( JT) 

RDYGN ( 1 ) =RDYGN ( JB ) 

RDYGN (NYP2) =RDYGN ( JT) 

C 

DO  590  J=1 7  NYP1 

DYGNH ( J) =0. 5* (DYGN ( J) +DYGN (J+l ) ) 

590  RDYGNH (J) =0. 5* (RDYGN (J) +RDYGN (J+l ) ) 

C 

RETURN 
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ENTRY  OGRID  <XGN*  YGN) 

RESET  OLD  GRID  PARAMETERS,  IN  PREPARATION  FGR  A  NEW  TIME  STEP 
IF (.NOT. XCHNG)  GO  TO  593 

DO  592  1=1 , NXP2 

592  DXGO  < I ) =DXGN ( I ) 

593  CONTINUE 

IF (.NOT. YCHNG)  GO  TO  595 

DO  594  J=1 , NYP2 

594  D YGO ( J ) =DYGN ( J ) 

GO  TO  595 

ENTRY  ORIGRD(XGN, YGN) 


ORIGINATE  THE  GRID 

SET  DEFAULT  :  GRID  IS  MOVING 
XCHNG=. TRUE. 

YCHNG=. TRUE. 

595  CONTINUE 

IF (.NOT. XCHNG)  GO  TO  597 

DO  596  1  =  1 »  NXF'l 

596  XGQ  < I ) =XGN ( I ) 

597  CONTINUE 

IF(. NOT. YCHNG)  RETURN 

DO  598  J=1 ,  NYF'l 

598  YGO ( J) =YGN ( J) 

RETURN 
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ENTRY  VOLFLX (U, V, DT) 

C  - - 

C 

C  EVALUATE  X  AND  Y  VOLUMETRIC  FLUX  THROUGH  INTERFACES 

C 

C 

DT2=0. 5*DT 
C 
C 

DO  602  J=i i NY 
DO  602  1  =  1 »  NXP1 

602  ADUDT ( I , J ) =U ( I , J+ 1 ) +U ( I + 1 , J+ 1 ) 

C 

DO  604  J=1,NYP1 
DO  604  1=1, NX 

604  ADVDT < I , J) =V  < 1  +  1 , J) +V  < 1  +  1 , J+l ) 

C 

GO  TO  (610,620,630,640,650,660,670)  GEOM 
C 
C 

610  CONTINUE 
C 

C  CARTESIAN  COORDINATES 

C 

DO  611  J=1 , NY 

DO  611  1=1 , NXP1 

ADUDT ( I , J) =ADUDT ( I , J) *DT2 

ADUDT ( I , J) =ADUDT < I , J) -DXG  ( I ) 

611  ADUDT  < I , J) =ADUDT ( I , J) *AY  < J) 

C 

DO  612  J=1 , NYP1 

DO  612  1=1, NX 

ADVDT ( I , J) =ADVDT ( I , J) *DT2 

ADVDT ( I , J) =ADVDT < I , J) -DYG  <  J ) 

612  ADVDT ( I , J) =ADVDT ( I , J) *AX ( I ) 

C 

RETURN 

C 

620  CONTINUE 
C 

C  CYLINDRICAL  R-Z  COORDINATES 

C 

DO  621  J=1 , NY 
DO  621  1  =  1 >  N  X P 1 

ADUDT ( I , J) =ADUDT ( I , J) *S0  < I ) *DT2 
ADUDT ( I , J) =ADUDT ( I , J) -DXG ( I ) 

621  ADUDT ( I , J) =ADUDT ( I , J) *AY(J> 

C 

DO  622  J= 1 , NYP1 

DO  622  1=1, NX 

ADVDT ( I , J) = ADVDT ( I , J) *DT2 

ADVDT  < I , J) =ADVDT ( I , J) -DYG ( J) 

622  ADVDT ( I , J) =ADVDT ( I , J) *AX (I ) 

C 

RETURN 

C 
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630  CONTINUE 

C  CYLINDRICAL  R-FYE  COORDINATES 

640  CONTINUE 

C  CYLINDRICAL  Z-FYE  COORDINATES 

650  CONTINUE 

C  SPHERICAL  R-THETA  COORDINATES 

660  CONTINUE 

C  SPHERICAL  R-FYE  COORDINATES 

670  CONTINUE 

SPHERICAL  THETA-FYE  COORDINATES 

PRINT  70 
STOP 


onn  n  n  non  nnn  n  nnn  non  n  n  n  n  nnnnn 


ENTRY  SORCES ( SRCTYP , SORCE , DT ) 


MANAGEMENT  OF  SOURCE  TERM  EVALUATION 


IF (SRCTYP. EQ.BDF)  GO  TO  750 

IF (SRCTYP.  EO. XGRD)  GO  TO  (760,300,999,999,999,999,999)  GEOM 
IF (SRCTYP. EQ. YGRD)  GO  TO  (780,780,999,999,999,999,999)  GEOM 


IF (SRCTYP. EQ. DIV)  GO  TO  (760,950,999,999,999,999,999)  GEOM 

TEXT (6) =TEXT (18) 

TEXT (7) =TEXT (19) 

TEXT (8) =TEXT (20) 

PRINT  10,  (TEXT(I) , 1=1, 11) 


STOP 


ENTRY  CLRSRC 


CLEAR  SOURCE  TERM 
LSRC=. FALSE. 

DO  710  J= 1 , NY 
DO  710  1=1, NX 
710  SOURCE ( I, J)=0. 

RETURN 


ENTRY  BODY (SORCE, DT) 


EVALUATE  BODY  FORCE  TYPE  SOURCE  TERMS 
750  CONTINUE 

LSRC=. TRUE. 

DO  755  J= 1 , NY 
DO  755  1=1, NX 
TEMPS (I, J)=AX (I) *AY(J) 

TEMP4 ( I , J) =SORCE ( 1+1 , J+l ) #DT 
TEMF‘5  ( I ,  J)  =TEMP3  ( I ,  J)  *TEMP4  ( I ,  J) 
755  SOURCE ( I , J) =SOURCE ( I , J) +TEMP5 ( I , J) 

RETURN 


onnn  n  n  nnnnnnnn  n  n  nnn 


T  ^ 


ENTRY  XGRAD ( SQRCE  >  DT ) 


EVALUATE  CARTESIAN  X  GRADIENT  COMPONENT 
SRCTYP=XGRD 
760  CONTINUE 

LSRC=. TRUE. 

DT2=0. 5#DT 

DO  765  J= 1 i NY 
DO  765  1  =  1 >  NXP1 

765  TEMPS ( I ,  J) =SORCE < I »  J+l ) +SORCE ( 1+1 , J+l ) 

DO  770  J=1 »  NY 
DO  770  1=1 i NX 

TEMP4 ( I , J) =TEMP3 ( 1+1 , J) -TEMPS C I ,  J) 
TEMPS ( I ,  J) -TEMP4 <  I ,  J) *AY ( J) *DT2 
770  SOURCE ( I , J ) =SQURCE (I ,  J ) +TEMP5 ( I ,  J) 

IF (SRCTYP. EQ. DIV)  GO  TO  780 

RETURN 


ENTRY  YGRAD ( SORCE , DT ) 


EVALUATE  CARTESIAN  Y  GRADIENT  COMPONENT 
790  CONTINUE 

LSRC=. TRUE. 

DT2=0.5*DT 


DO  785  J= 1 i NYP 1 
DO  785  1=1, NX 

785  TEMPS ( I , J) =SORCE  < 1+1 , J) +SORCE ( 1  +  1 , J+l ) 

DO  790  J=1 , NY 
DO  790  1=1, NX 

TEMP4 ( I , J) =TEMP3 ( I , J+l ) -TEMPS ( I , J) 
TEMPS  < I , J) =TEMP4 ( I , J) *AX  < I ) *DT2 
790  SOURCE  < I , J) =SOURCE ( I , J) +TEMP5 ( I , J) 

RETURN 
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ENTRY  RCGRAD ( SORCE  ,  DT ) 


C  - 

C 

C  EVALUATE  CYLINDRICAL  R  GRADIENT  COMPONENT 

800  CONTINUE 

LSRC=. TRUE. 

C 

DO  805  J=1,NY 
DO  805  1  =  1 ,  NXP1 

TEMPS ( I , J ) =SORCE ( I , J+ 1 ) + SORCE ( I + 1 , J+ 1 > 
TEMF'4  ( I  ,  J)  =SG  ( I )  #  AY  <  J )  *DT 
805  TEMPS (I, J>=0.5*TEMP3(I, J) 

C 

DO  810  J=1 ,  NY 
DO  810  1=1, NX 

TEMPS ( I , J) =TEMP3 < I , J) *TEMP4 ( 1*1 , J) 

810  TEMPS ( I ,  J) =TEMP4 ( 1+1 , J) -TEMP4 ( I , J) 

C 

DO  815  J= 1 , NY 
DO  815  1=1, NX 

TEMP4  ( I ,  J)  =TEMF'3  ( 1  + 1 ,  J )  -TEMPS  ( I ,  J) 
TEMPS ( I , J) =TEMP5 ( I , J) *SORCE (1+1 , J+l ) 
SOURCE ( I , J) =SOURCE ( I , J) +TEMP4 ( I , J) 

815  SOURCE ( I , j) =SOURCE ( I , J) +TEMP5 ( I , J) 

C 


RETURN 

C 

C 

C 

ENTRY  RCDIV (SORCE, DT) 

C  - 

c 

C  EVALUATE  CYLINDRICAL  DIVERGENCE 

950  CONTINUE 

LSRC=. TRUE. 

DT2=0. 5*DT 
C 

DO  955  J= 1 , NY 
DO  955  1  =  1 ,  NXF'l 

TEMPS ( I , J ) =SORCE ( I , J+ 1 ) +SORCE ( I + 1 , J+ 1 ) 
955  TEMP4 < I , J) =SQ ( I) *AY ( J) *0T2 


960 

C 

C 

C 

999 

C 

C 

C 


DO  960  J=1 , NY 
DO  960  1=1, NX 

TEMPS  ( I,  J)  =TEMF'S  ( I  + 1 ,  J)  -TEMPS  ( I ,  J) 
TEMPS  ( I ,  J)  =TEMP5  ( I ,  J)  *  TEMF'4  ( 1  +  1  ,  J) 
SOURCE ( I , J) =SOURCE ( I , J ) +TEMP3 < I , J) 

GO  TO  780 

RETURN 

PRINT  70 
STOP 
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ENTRY  SETGOM <  GOMTRY , CRD 1 , CRD2 , N 1 , N2  > 


SET  AND  CHECK  REQUEST  FOR  A  PARTICULAR  GEOMETRY 


GE0M=0 

IF (GOMTRY. ED- TGM < 1) )  GE0M=1 

IF (GOMTRY. NE. TGM (2) )  GO  TO  1210 
IF  (CR'Dl  .NE.  TCRD  (4)  )  GO  TO  1205 
IF (CRD2. EO. TCRD (3) )  GE0M^2 
IF (CRD2. EQ. TCRD (6) )  GE0M=3 
1205  CONTINUE 

IF (CRD1.EQ. TCRD ( 3 ) . AND. CRD2. EQ. TCRD (6) )  GE0M=4 

1210  CONTINUE 

IF (GOMTRY. NE. TGM (3) )  GO  TO  1220 
IF (CRD1 . NE. TCRD (4) )  GO  TO  1215 
IF (CRD2. EQ. TCRD (5) )  GE0M=5 
IF (CRD2. EQ. TCRD (6) )  GEQM=6 
1215  CONTINUE 

IF (CRD1 . EQ. TCRD (5) . AND. CRD2. EQ. TCRD (6) )  GEGM=7 

1220  CONTINUE 

IF(GEOM.GT.O)  GO  TQ  1225 

ISSUE  AN  ERROR  MESSAGE  UPON  REQUEST  OF  AN  UNRECOGNIZED  GEOMETRY 
AND  STOP 

TEXT (6) =TEXT (12) 

TEXT  <7) =TEXT  < 13) 

TEXT (8) =TEXT ( 14) 

PRINT  10,  <TEXT(I) , 1=1, 11) 

STOP 

1225  CONTINUE 


C 


C 


C 


C 

C 


NX=N  1 
NXP 1 =NX+ 1 
NXF'2=NX+2 


NY=N2 
NYP 1 =NY+ 1 
NYP2=NY+2 

XPRDC=. FALSE. 

IL=2 

IR=NXF'l 

YF'RDC  =  .  FALSE. 

JB=2 
JT =NYP 1 

LSRC=. FALSE. 
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C  SET  DEFAULT  :  ALL  BOUNDARIES  ASSUMED  PERMEABLE 

DO  1226  I  — 1  *  NXP2 
PRMBLB ( I ) =1 . 0 
1226  PRM8LT  < I ) — 1 . 0 

C 

DO  1227  J=1  >  NYF’2 
PRMBLL  < J) =1 . 0 


1227  PRMBLR ( J) =1 . 0 

C 


nnnn  nn  n  nn  nnonnn 


ENTRY  PROD I C ( CRDNT  »  CRD ) 


IDENTIFY  X  OR  Y  COORDINATE  AS  PERIODIC 


IF (CRDNT. NE . 1 )  GO  TO  1230 

XPRDC=. TRUE. 

IR=2 

1L=NXF1 

PRINT  40,  CRD 

1230  CONTINUE 

IF (CRDNT. NE. 2)  GO  TO  1240 

YPRDC=. TRUE. 

JT =2 
JP=NYP1 

PRINT  40,  CRD 

1240  CONTINUE 

IF (CRD.  EQ. TCRD (4) . OR. CRD. EO. TCRD  (5) )  PRINT  20, CRD 
RETURN 
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1251 


1252 


1255 


1260 


/ 


ENTRY  SETLMT (LMTR1 ,  LMTR2) 


SET  AND  CHECK  REQUEST  FOR  A  PARTICULAR  FLUX  LIMITER 


DO  1255  1=1,4 
DO  1251  J=1 , 2 

IF(LMTR1 (J) . NE. TLM1 < J, I) >  GO  TO  1255 
CONTINUE 

DO  1252  J= 1 , 2 

IF (LMTR2 ( J) . NE. TLM2  < J, I ) )  GO  TO  1255 
CONTINUE 

GO  TO  1260 

CONTINUE 

ISSUE  AN  ERROR  MESSAGE  UPON  REQUEST  OF  AN  UNRECOGNIZED  FLUX 
LIMITER  AND  STOP 

TEXT (6) =TEXT (15) 

TEXT (7) =TEXT (16) 

TEXT (8) =TEXT  < 17) 

PRINT  10,  <TEXT(I) ,  1  =  1, 1  1) 

STOP 

CONTINUE 

ILMTR=I 

RETURN 
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ENTRY  SOLDFY (BONDRY, KSTRT, KEND) 


CANCEL  THE  DIFFUSION  AND  ANTI-DIFFUSION  FLUXES  THROUGH  F'ATCHE 
OF  BOUNDARY  INTERFACES 


IF (BONDRY. NE.TBND(l) )  GO  TO  1280 

DO  1275  J=KSTRT , KEND 
PRMBLL ( J) =0. 

PRINT  60,  BONDRY, KSTRT, KEND 
CONTINUE 

IF (BONDRY. NE. TBND (2) )  GO  TO  1290 

DO  1285  J=KSTRT , KEND 
PRMBLR ( J) =0. 

PRINT  60,  BONDRY, KSTRT, KEND 
CONTINUE 

IF (BONDRY. NE. TBND (3) )  GO  TO  1300 

DO  1295  I=KSTRT , KEND 
PRMBLB ( I ) =0. 

PRINT  60,  BONDRY, KSTRT, KEND 
CONTINUE 

IF (BONDRY. NE. TBND (4) )  GO  TO  1310 

DO  1305  I*=KSTRT , KEND 
PRMBLT ( I ) =0. 

PRINT  60,  BONDRY, KSTRT, KEND 


CONTINUE 

RETURN 


non 


SUBROUTINE  NUMU  <NI , N J i EPS »  NUV >  MUV ) 


EVALUATE  DIFFUSION  AND  ANTI-DIFFUSION  COEFFICIENTS 


C 


PARAMETER- 

PARAMETER 

PARAMETER 


NPX=100,NPY=100 
NP1X=NPX+1 ,NP1Y=NPY+1 
NF‘2X=NPX+2,  NP2Y=NPY+2 


REAL  EPS  <  NP 1 X  ,  NP 1 Y ) , NUV <NP2X , NP2Y) ,MUV(NP2X 


C 

DO  100  J= 1 i NJ 
DO  100  1=1, NI 

EPS ( I , J) =EPS ( I , J) *EPS <  I ,  J) 
NUV ( I , J) =0. 333333*EPS ( I , J) 
NUV ( I , J) =0. 166667+NUV ( I , J) 
100  MUV ( I , J ) =NUV ( I , J ) -EPS ( I , J ) 

C 

RETURN 

END 
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